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Abstract 

The  objective  of  this  paper  is  to  establish  a  simple  two- 
dimensional  theoretical  model  in  an  attempt  to  use  a  computer  to 
mwerically  simulate  the  experimental  results  of  Hummel  regarding 
the  rolling-up  and  interaction  of  the  leading-edge  and  trailing- 
edge  vortex  sheets  on  a  delta  wing.  It  was  found  experimentally 
that' when' the  leading  vortex  is  present  the  trailing-edge  vortex 
sheet  will  roll  up  another  vortex  downstream  from  the  trailing- 
edge.  Furthermore,  the  circulation  of  the  leading-edge  vortex  is 
opposite  in  direction  to  that  of  the  trailing-edge  vortex.  The 
numerical  results  are  in  good  agreement  with  the  experimental 
pictures. 

I.  Introduction 

One  of  the  problems  of  major  concern  for  researchers  in 
aerodynamics  and  aircraft  designers  is  the  non-linear  aerodynamic 
characteristics  caused  by  the  separation  of  body  and  wing 
vortices  of  the  aircraft  at  large  angles  of  attack.  Effective 
utilization  of  the  additional  lift  generated  by  body  and  wing 
vortices  can  improve  the  aerodynamic  properties  of  the  aircraft 
and  increase  the  maneuverability.  To  study  the  mechanism  of 
formation  of  body  and  wing  vortices  and  the  complicated 


interaction  between  various  vortex  systems  as  well  as  between 
vortices  and  the  aircraft  in  detail  is  the  key  to  the  accurate 
estimation  of  various  non-linear  force  and  torque  terms  on  the 
aircraft.  Therefore,  the  study  of  vortex  motion  has  important 
practical  values. 

The  study  of  the  leading-edge  vortex  of  a  slender  delta  wing 
began  in  the  forties  and  fifties.  There  are  significant  advances 
in  recent  years.  In  addition  to  measuring  force  and  pressure, 
recent  experimental  studies  focused  on  the  application  of  display 
technique  to  the  flow  field  as  well  as  on  the  detection  of  fine 
details  of  the  spatial  flow  field.  Based  on  the  "contours"  of 
total'  pressure,  static  pressure  and  dynamic  pressure  measured,  as 
well  as  on  the  spatial  distribution  of  the  flow  direction,  we  can 
have  a  more  direct  and  profound  understanding  of  the  vortex  flow 
field. 

In  the  early  stage,  the  theoretical  study  of  leading-edge 

vortex  was  based  on  the  conic  flow  assumption  which  simplified  a 

three-dimensional  flow  problem  to  a  two-dimensional  problem  on  a 

transverse  plane,  including  the  work  done  by  C.E.  Brown  and  W.H. 

Michael^ K.W.  Mangier  and  J.H.B.  Smlth^^,  and  the  later 

r  3 1 

improvement  made  by  J.H;B.  Smith  .  The  Smith  model  divides  the 
vortex  layer  into  two  points.  The  outer  part  uses  a  broken  line 
section  to  replace  the  vortex  layer.  The  inner  part  uses  a 
concentrated  vortex  to  represent  the  core  and  a  "vortex  transport 
line"  to  connect  inner  and  outer  regions.  The  model  could  be 
used  to  obtain  the  shape  of  the  vortex  layer,  and  the  strength 
and  intensity  of  the  core.  However,  the  accuracy  is  not 


desirable.  Alter  computers  are  extensively  used,  a  vortex 
lmt&iee  — thod  with  leading-edge  separation  vortex,  introduced  by 
C.M.BetouepkoBckHH^^  and  O.A.  Kandil,  D.T.  Mook  and  A.H. 
Wayfeh^^ ,  la  a  representative  aethod.  The  leading-edge  vortex 
layer  la  replaced  by  several  discrete  vortex  threads.  Through 
iterations,  the  position  of  free  vortex  threads  are  determined. 
The  boundary  conditions  on  the  wing  surface  are  also 
simultaneously  satisfied.  In  order  to  accurately  calculate  the 
load  distribution  on  the  wing,  P.E.  Rubbert  et  al^^  introduced 
the  "free  vortex  layer"  method  by  using  higher  order  surface 
elements.  It  can  be  used  to  calculate  the  shape  of  non-conical 
flow- fields  and  vortex  layers,  as  well  as  the  load  distribution 
on  the  entire  wing. 

All  the  experimental  and  theoretical  studies  discussed  above 
are  focused  on  the  rolling-up  of  the  leading-edge  vortex  layer, 
the  force  and  torque  characteristics,  and  the  calculation  of  load 
distribution.  It  seems  that  there  is  little  work  done  on  the 
development  of  leading-edge  vortex  at  downstream  from  the 
trailing-edge  and  the  interaction  between  leading-edge  and 
trailing-edge  vortices.  We  are  very  much 
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intrigued  by  the  work  done  by  D.  Hummel^®^.  Hummel  performed  a  ^ 
series  of  fine  manuscripts.  In  particular,  he  did  an 
experimental  study  of  the  interaction  between  leading-edge  and 
trailing-edge  vortices.  From  his  measured  total  pressure,  static 
pressure  and  spatial  flow  direction  distribution,  we  can  see  that 


two  spiral  vortices  are  gradually  formed  downstream  from  the 
tr ailing-edge .  One  is  the  leading-edge  vortex  and  the  other  Is 
the  vortex  rolled-up  by  the  trailing  vortex  layer.  The 
circulations  of  these  two  vortices  are  opposite  in  direction.  A 
schematic  diagram  of  the  flow  pattern  is  shown  in  Figure  1. 


Figure  1.  Schematic  Diagram  for  the  Formation  of  Downstream 
Vortices  of  a  Slender  Delta  Wing 

Inspired  by  Hummel's  experiment  results,  we  attempted  to 
establish  a  simple  theoretical  model  to  simulate  Hummel's  results 
numerically  on  a  computer.  This  study  will  benefit  the 
understanding  of  the  structure  of  a  down  wash  flow  field. 

II.  Theoretical  Analysis 

In  order  to  study  the  interaction  between  leading-edge  and 
trailing-edge  vortices,  we  must  first  obtain  the  shape,  position 
and  strength  of  the  rolling-up  of  the  leading-edge  at  the 
trailing-edge.  In  addition,  we  must  also  have  the  intensity 
distribution  of  the  trailing-edge  vortex,  i.e.,  the  vortex 


intensity,  or  spanwise  circulation,  distribution  on  the  wing. 
Figures  2  and  3  show  the  pressure  distribution  on  the  wing 
surface  and  the  vortex  line  shape  measured  by  Hummel.  From  the 
figures,  one  can  see  that  the  surface  pressure  distribution  and 
the  vortex  line  are  essentially  different  from  those  obtained 
based  On  the  linearized  slender  wing  theory  of  Jones  due  to  the 
presence  of  the  leading-edge  vortex.  However,  as  compared  to 
Smith's  theory,  the  shape  of  the  pressure  distribution,  the 
position  of  the  suction  peak  and  the  shape  of  the  vortex  are 
qualitatively  similar.  However,  there  are  some  differences 
quantitatively.  In  other  words,  as  a  preliminary  theoretical 
investigation,  a  two-dimensional  model  can  reflect  the  major 
characteristics  of  the  flow  field.  But,  we  did  not  choose 
Smith's  vortex  layer  model.  Instead,  a  simpler  two-dimensional 
unsteady  flow  analogy  was  used.  Our  theoretical  model  was  built 
based  on  a  discrete  vortex  method,  which  does  not  require 
iterations  to  solve  a  set  of  non-linear  equations. 


Figure  3.  Adhered  vortex 
vector  (right)  and  adhered 
vortex  line  (left) 


Figure  2.  Pressure  distri¬ 
bution  on  Delta  wing  A»1 
o*20 .5 

1 — experimental 


Figure  4 
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Coordinate  systems  definition 


It  is  an  initial  value  problem  of  a  series  of  normal  differential 
equations.  It  is  not  limited  to  the  "conic  flow"  assumption 
which  will  facilitate  the  extension  to  more  complicated  airfoils. 
It  also  facilitates  the  further  consideration  of  "secondary 
vortex"  separation  problems. 

Based  on  the  two-dimensional  unsteady  flow  analogy,  the 
three-dimensional  flow  of  a  delta  wing  with  an  attack  angle  can 
be  considered  as  an  unsteady  flow  around  a  two-dimensional  plate 
in  the  x  plane.  The  width  of  the  plate  at  any  time  corresponds 
to  the  wing  span  in  the  x  position.  When  the  flow  passes  the 
edge  of  the  plate,  the  boundary  layer  is  separated.  A  free  shear 
layer  is  -formed  due  to  the  velocity  difference  between  the  upper 
and  lower  surface.  Based  on  existing  studies,  it  is  known  that, 
as  long  as  it  does  not  involve  the  mechanism  of  separation  (which 
is  a  viscosity  effect),  the  free  shear  layer  after  separation  can 
be  assumed  as  an  inviscid  vortex  layer.  As  a  next  step,  the 
vortex  layer  is  replaced  by  several  discrete  point  vortices. 
Therefore,  the  flow  around  the  plate  satisfies  LaPlace  equation 

~  ~  ( 1 ) 
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From  this  point  on,  the  plate  can  be  transformed  into  a  circle  by 
using  a  complex  function  method.  Hence,  the  mathematics  of  the 
problem  becomes  a  flow  around  a  cylinder  with  a  finite  number  of 
point  vortices  outside  the  circle.  In  this  case,  the  complex 
potential  expression  of  the  flow  is: 


where  the  first  term  is  the  complex  potential  of  a  uniform 
incident  flow  around  a  cylinder  and  the  second  term  is  the 
complex  potential  generated  by  the  finite  number  of  point 
vortices  and  their  image  vortices  outside  the  circle.  The 
boundary  conditions  are  automatically  satisfied  on  the  circle. 
The  conformal  mapping  is  s 


X  - 


(3a) 


f  -  X  +  Vx  -  /» 


(3b) 


It  should  be  noted  that  in  this  transformation  there  is  a 
magnifying  factor  dX/dtj^^  at  infinity.  Therefore,  IT  =  l/2V°°sina. 
The  diameter  of  the  circle  is  the  width  of  the  plate. 

In  addition  to  surface  boundary  conditions,  the  Kutta 
condition  must  also  be  satisfied. 


Because  dC/dX 
dW(?)/d?  *0. 
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••  at  the  edge  of  the  wing,  therefore,  we  must  have 
Hence,  we- get 
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The  point  vortex  moves  downstream  at  the  local  velocity. 
Therefore,  we  must  also  find  the  velocity  of  every  point  vortex. 


If  the  velocity  of  the  kth  vortex  is  expressed  as: 


4W(X) 


*-1,2, 


where  the  second  term  on  the  right  indicates  that  the  inducing  /457 
velocity  of  the  kth  vortex  itself  should  be  subtracted  from  the 
calculation  of  the  velocity  at  the  kth  vortex  because  it  is  a 
velocity  singular  point.  After  using  equations  (2)  and  (3)  to 
calculate,  equation  (6)  becomes 


where 


ft  —  —  (G,  +  Gj  +  GO  |  +  G« 
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Intensity  and  Position  of  Newly  Generated  Leading-edge 
Vortex 

In  the  unsteady  flow  analogy,  the  leading-edge  vortex  is 
approximated  by  many  discrete  vortices.  Newly  generated  vortices 
continue  to  be  separated  from  the  edge  of  the  wing  into  the  flow 
field  with  time.  Therefore,  the  number  of  vortices  continues  to 
increase  in  the  flow  field.  The  intensity  of  a  newly  generated 
vortex  has  a  great  effect  on  the  shape  and  position  of  the 
rolling-up  of  the  leading-edge  vortex  and  the  surface  pressure 
distribution.  Many  authors  have  investigated  this  problem.  The 


variation  of  the  vortex  flow  in  the  shear  layer  near  a  leading- 
edge  separation  point  with  time  is: 

<"•/«-  ^down-V  <8) 

where  Vup  and  V^own  represent  the  velocities  at  the  upper  and 
lower  surface  of  the  shear  layer  near  a  separation  point. 
According  to  Sacks^^  method,  let  Vg=%(Vup  an<*  Vdown^*  !•*•»  the 
average  velocity  of  the  shear  layer  which  is  also  the  velocity  at 
which  the  shear  layer  is  dragged  out  the  edge  of  the  wing,  y 
-Vup  is  the  vortex  intensity  on  a  unit  length.  Based  on  these, 
equation  (8)  can  be  re-written  as: 

Ar  “At  .  (Vgyx)  -  (VgAt)  .  yx»AS  .  yx  (9) 

As  is  the  length  of  the  shear  length  dragged  out  in  At  time.  The 
expression  for  Vg  is  found  to  be 

■  I  y  -r.  [  t! _ E_1  <10) 
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Because  of  symmetry,  w  is  equal  to  zero  in  practice.  There  is 
only  a  y-direction  velocity  which  shows  that  the  shear  layer  is 
dragged  out  tangentially  from  the  leading-edge.  In  our 
computation,  v  was  found  with  equation  (10)  in  order  to 

9 

determine  the  vortex  layer  length  As.  Then,  the  simplified  point 
vortex  is  centered  in  this  section  of  vortex  layer  so  that  the 
intensity  of  the  newly  generated  vortex  could  be  determined  by 
using  the  Kutta  condition. 

Surface  Pressure  Distribution  and  Attached  Vortex  Line 
From  the  definition  of  the  pressure  coefficient  c  =P-P  / 


*  -*»*•  -  jf*  eow- 


(11) 


where  9  »v,  •  «w ,  f»R.P.[W].  The  calculation  of  »„  must  take  two 

y  *  * 

aspects  into  account.  One  is  that  the  semi-wing  span  s(x)  is  a  / 4 
functionof  x.  The  variation  caused  by  s(x)  is  »x»ds(x)/dx.R.P.0 
W(X)/as .  The  other  is  *x  caused  by  the  variation  of  the  point 
vortices  with  x.  The  derivation  of  the  entire  formula  is 
tedious.  It  is  omitted  here.  On  the  wing  surface,  *x-0.  The 
attached  vortex  line  on  the  wing  surface  can  be  determined  based 


on  the  following 


*-n  x  (VUp“Vdown) 


(12) 


Y  is  the  attached  vortex  vector  on  the  wing  surface:  Let  yx  and  Ty 
be  the  components  of  t  in  the  x  and  y  directions,  respectively. 

^x^+^y^  ‘  Then,  *  Yy“*xup”*xdown’  w^ere  *xup’ 

f xdown ’  fyup»  and  fydown  are  the  velocity  components  on  the 
upper  and  lower  wing  surface. 

After  the  attached  vortex  line  on  the  wing  surface  is  found, 
it  is  very  easy  to  obtain  the  intensity  of  the  trailing-edge  tail 
vortex.  Due  to  the  fact  that  *rxs-0r/8y,  the  intensity  of  the 
i**1  tail  vortex  as  a  result  of  trailing-edge  discretization  is 

^ Ari^trailing-edge-/  i+^irx^ trail ing-edgedy 

yi 

III.  Brief  Description  of  Results  of  Computation 
Mathematically,  this  computation  is  to  solve  a  variable 
number  of  first  order  normal  differential  equations.  We  can  use 
the  Runge-Kutta  method  to  calculate  gradually  from  the  apex  of 


WvVv';-.'.. 


the  wing  downstream.  With  each  increasing  step,  a  new  vortex  is 
generated , 


Figure  6.  Surface  Pressure  Distribution  and  Attacned  Vortex 
Distribution  (Trailing-edge  of  the  Wing) 


resulting  in  two  additional  equations.  When  calculating  the  / 459 

trailing-edge ,  the  intensity  of  the  trailing-edge  vortex  is  found 
and  combined  into  the  original  equation  to  be  moved  downstream. 

In  order  to  compare  with  the  experiment,  we  calculated  a 
delta  wing  whose  aspect  ratio  A»1.0  and  attack  angle  a=20.5°. 

Figure  5  shows  the  gradual  rolling-up  of  the  leading-edge  vortex 
over  the  upper  wing  surface.  Figure  6  shows  the  pressure 
distribution  and  attached  vortex  intensity  distribution  at  the 
trailing-edge.  The  pressure  distribution  is  very  close  to  that 
calculated  by  Smith.  However,  it  is  different  from  the 
experimental  result  (see  Figure  2).  From  the  experimental  result 
one -can  see  that  the  flow  in  the  front  part  of  the  delta  wing 
approaches  the  conical  flow  assumption.  However,  the  rear  part, 
especially  near  the  trailing-edge,  is  no  longer  a  conical  flow. 

The  suction  peak  decreases  with  Increasing  x.  But,  this  tendency 
cannot  be  calculated  using  a  two-dimensional  model.  This  is 
because  the  two-dimensional  model  does  not  meet  the  trailing-edge 
Kutta  condition.  Although  the  load  distribution  on  the  wing 
surface  can  be  more  accurately  calculated  based  on  a  three- 
dimensional  flow  model  using  a  higher  order  surface  element  "free 
vortex  layer"  method  currently  under  development,  yet  it  takes 
too  much  computing  time.  As  a  qualitative  analysis,  we  chose  the 
two-dimensional  model. 

From  the  distribution  of  attached  vortex  intensity  along  the 
span  yx  one  can  see  that  y  is  negative  over  most  of  the  wing 
span.  It  is  positive  near  the  edge  of  the  wing.  In  addition, 
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under  the  leading-edge  vortex,  y%  has  a  negative  maximum. 
According  to  the  calculated  result  reported  in  reference  [9],  a 
vortex  will  be  rolled-up  at  the  extremum  |yxI*  When  yx  is 
negative,  the  vortex  rolled-up  is  clockwise. 

Figure  7  shows  the  rolling-up  of  leading-edge  and  trailing- 
edge  vortices  and  their  interaction.  For  comparison,  Hummel's 
experimental  results  are  again  shown  in  Figure  8.  From  the 
figures,  the  two  situations  are  quite  similar.  In  Figure  7(a), 
the  tralllng-edge  vortex  layer  already  begins  to  fluctuate.  It 
bulges  slightly  at  the  extremum  |yx|and  develops  downstream.  On 
one  hand,  it  continues  to  bulge  and  enlarge  and  gradually  rolls 
up  into  &  clockwise  vortex.  On  the  other  hand,  because  of  the 
side  wash  velocity  effect  induced  by  the  leading-edge  vortex,  the 
traillng-edge  vortex  layer  extends  in  the  direction  of  the  wing 
edge.  The  vortex  rolled  up  by  the  trailing-edge  moves  outward. 

It  is  initially  on  the  right  lower  side  of  the  leading-edge 
vortex  and  then  gradually  rises.  From  the  figure  one  can  also 
see  that  the  trailing-edge  vortex  begins  to  roll  up  at 
approximately  1/4  of  a  wing  span  (x/cr«1.10)  from  the  trailing- 
edge.  At  1/2  wing  span  (x/cr*1.20),  it  has  already 
developed  well. 


Figure  7.  Rolling-up  of  Leading-edge  and  Trailing-edge  Vortex 
Layers  Downstream  from  Trailing-edge  of  the  Wing 


In  order  to  study  the  effect  of  the  leading-edge  vortex  on 
the  trailing-edge  vortex,  we  also  did  two  Interesting  numerical 
experiments.  Figure  9  shows  the  effect  of  the  leading-edge 
vortex  layer.  We  artificially  neglected  the  vortex  layer  and 
consolidated  the  leading-edge  as  a  point  vortex.  The 
consolidation  is  based  on  the  conservation  of  vortex  moment  and 
circulation.  In  the  figure,  the  symbol  a  represents  the 
consolidated  leading-edge  vortex.  We  found  in  the  figure  that 
the  trailing-edge  vortex  could  also  roll-up  a  vortex.  However, 
the  shape  and  position  are  quite  different  from  those  shown  in 
Figures  7  and  8.  Thus,  the  effect  of  the  leading-edge  vortex 
layer  cannot  be  neglected. 


Figure  9.  Rolling-up  of  Trailing-edge  Vortex  When  Leading-edge 
Vortex  is  Consolidated  into  a  Point  Vortex 
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Figure  10.  Rolling-up  of  Trailing-edge  Vortex  When  Neglecting 
Leading-edge  Vortex 

Figure  10  shows  the  rolling-up  of  the  trailing-edge  vortex  /461 
when  the  entire  leading-edge  is  neglected.  Just  as  expected,  a 
clockwise  vortex  is  rolled-up  at  the  extremum  |rx|.  A 
counterclockwise  vortex  is  rolled-up  at  the  wing  tip.  Because 
the  rx  value  is  very  small  at  the  wing  tip,  only  a  small  vortex 
is  rolled-up.  In  addition,  the  number  of  point  vortices  is  not 
sufficient  to  see  clearly.  Because  of  the  absence  of  the 
leading-edge  vortex  side  wash  velocity  effect,  the  trailing-edge 
vortex  layer  extends  slowly.  The  clockwise  vortex  is  on  the 
Inside  of  the  wing  tip  vortex. 


IV.  Conclusions 


In  this  work,  a  two-dimensional  discrete  vortex  model  was 
established  based  on  the  two-dimensional  unsteady  analogy. 
Numerical  simulation  of  Hummel's  experimental  results  was 
realized  on  a  computer.  The  rolling-up  of  the  leading-edge  and 
trailing-edge  vortices  and  the  results  of  their  interaction  thus 
obtained  are  very  similar  to  Hummel's  experimental  results.  It 
proves  that  it  is  basically  feasible  to  study  the  mechanism  using 
a  two-dimensional  model.  Major  physical  pictures  of  the  flow 
field  can  be  obtained. 

1.  In  addition  to  the  vortex  rolled  up  by  the  leading-edge 
vortex  downstream,  the  trailing-edge  vortex  will  roll  up  another 
vortex.  The  circulations  of  these  two  vortices  are  opposite. 

2.  Under  the  influence  of  the  side  wash  velocity  induced  by 
the  leading-edge  vortex,  the  trailing-edge  vortex  layer  extends 
toward  the  edge  of  the  wing.  Initially,  a  vortex  is  rolled  up  on 
the  lower  right  side  of  the  leading-edge  vortex^  Then,  as  the 
circulation  gradually  increases,  it  rises  comparatively.  The 
presence  of  the  leading-edge  vortex  accelerates  the  rolling-up 
process  of  the  trailing-edge  vortex  and  also  pulls  it  outward. 

3.  The  wash  flow  field  is  complicated  where  there  are 
leading-edge  and  trailing-edge  vortices  present.  There  is  a  need 
to  understand  it  better.  This  study  has  helped  the  understanding 
of  the  physical  picture  of  the  wash  flow  field.  However,  because 
of  the  characterisic  deficiencies  of  the  two-dimensional  model, 
there  are  discrepancies  in  the  quantitative  determination  of  the 


pressure  distribution  on  the  wing  surface.  A  more  complex  three- 
dimensional  vortex  layer  model  must  be  used  to  more  accurately 
calculate  the  pressure  distribution  on  the  wing  surface.  This 
work  is  just  a  preliminary  investigation. 

After  this  paper  was  sent  for  review,  we  discovered  two 
similar  studies  done  abroad.  Kandil^^  used  a  vortex  lattice 
method  to  calculate  a  three-dimensional  flow  field.  But,  the 
structure  of  the  rolled-up  vortex  layer  is  not  clear.  The  work 
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done  by  Hoeijmakers  ’  et  al.  is  similar  to  ours.  They  also 
used  a  two-dimensional  vortex  layer  model  to  obtain  similar 
results.  However,  the  methodology  is  not  quite  the  same.  On  the 
wing'  surface  and  in  the  vortex  layer,  they  used  dipole 
distribution,  vortex  layer  shape  and  wing  surface  dipole  strength 
distribution  and  solved  them  by  iteration.  We  established  a 
series  of  point  vortex  equations  through  conformal  mapping  to 
convert  it  to  a  problem  of  solving  for  initial  values. 
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INVESTIGATION  OF  ROLLING-UP  AND  INTERACTION  OF 
LEADING-EDGE  AND  TRAILING-EDGE  VORTEX 
SHEETS  ON  A  SLENDER  DELTA  WING 


Tin  Xieyuan,  Xu  Nan,  Deng  Guohua 
(Umvtratj  of  Sdemcr  end  Technology  of  Chime) 


Abstract 


Hummel’s  experiment  on  the  rolling-up  and  interaction  of  the  leading-edge  and 
trailing-edge  vortex  sheets  at  slender  delta  wing  is  modeled  numerically  by  a  simple  two- 
dimensional  theory.  The  numerical  results  show  that  the  trailing-edge  vortex  sheet  will 
roll-up  at  the  downstream  of  the  wing,  in  the  presence  of  the  leading-edge  vortex,  and 
the  direction  of  the  rircuktion'of  the  leading-edge  rortex  is  oppsite  to  the  trialing- 
edge  vortex.  The  numerical  results  are  in  good  agreement  with  the  experiment.  This 
work  is  important  to  further  understand  of  the  downstream  flow-field  of  a  wing. 
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Calculation  of  Circular  Jet  with  Particles  Impacting  /463 

Upon  a  Plate 

Liu  Dayou 

(Institute  of  Mechanics,  Academia  Sinica) 

Abstract 

Under  the  inviscid  and  incompressible  flow  condition,  the- 

flow  field  of  a  uniform  circular  jet  with  particles  impacting 

upon  an  infinite  plate  is  calculated.  In  addition,  two  drag 

2/3 

coefficient  formulas,  i.e.,  24/Re  and  24/Re  (1+Re  ),  are  used 

to  calculate  the  trajectories  of  spherical  particles  in  the 
flow  field.  Assuming  particles  are  uniformly  distributed  in  the  jet 
outlet,  the  impact  coefficient  P(St,fi*)  curve  (known  as  the 
collection  probability  in  the  study  of  samplers)  has  been 
obtained.  The  rationale  of  each  assumption  is  discussed.  The 
effect  of  viscosity  is  discussed.  The  P(St,fi*)  curve  is 
corrected  for  the  effect  of  viscosity.  Results  indicate  that 
impact  points  are  mainly  concentrated  in  the  X<2  region  on  the 
plate.  Although  it  is  assumed  to  be  very  large  in  the 
calculation,  however,  the  results  are  still  applicable  to  the 
situation  H=1.5.  Experimentally,  it  has  been  proven  that 
samplers  designed  based  on  the  P(St,n*)  curve  calculated  in  this 
work  can  realize  anticipated  specifications. 

As  ecological  science  develops,  there  is  a  need  to  study  the 
effect  of  particles  of  various  diameters  in  the  atmosphere  on 
human  health.  Therefore,  -in  addition  to  the  need  to  know  the 
total  mass  density  and  number  density  of  particles  in  the 
atmosphere,  it  is  also  required  to  know  the  particle  diameter 


distribution.  To  this  end,  many  nations  are  developing  various 
atmospheric  particle  collectors  capable  of  sorting  by  diameter. 
One  of  the  most  common  types  is  the  impact  sampler  which  is  based 
on  fluid  dynamics  principles'^.  The  basis  for  studying  the 
collection  probability  of  an  impact  sampler  is  to  determine  the 
particle  trajectory  of  a  jet  impacting  upon  a  plate. 

Based  on  the  inviscid  and  incompressible  fluid  assumption, 
the  flow  field  of  a  uniform  jet  Impacting  upon  an  infinite  plate 
and  the  spherical  particle  trajectory  in  the  flow  field  are 
calculated.  The  impact  coefficients  P  of  various  diameters  (also 
know  as  collection  probabilities  in  the  study  of  samplers)  were 
determined.  The  rationale  for  each  assumption  made  in  the 
calculation  was  discussed.  Moreover,  some  corrections  were  made 
based  on  the  actual  flow  to  finally  obtain  the  P(St,n*)  curve  for 
the  design  of  samplers. 

I.  Basic-  Assumptions  and  Dimensional  Analysis 

If  a  gas  flow  is  injected  out  of  the  round  hole  CC  toward 
the  plate  AA  (as  shown  in  Figure  1),  the  streamline  (CB,PQ)  bends 
as  it  is  passing  through  the  plate.  Both  magnitude  and  direction 
of  velocity  change.  At* a  distance  from  the  axis,  such  as  the 
flow  along  the  plate  near  B  and  Q,  the  velocity  is  close  to  the 
exit  velocity  vm.  The  particles  in  the  flow  move  along  with  the 
gas  in  the  hole.  When  the  streamline  bends,  the  velocity  of  the 
particle  lags  behind  that  of  the  gas  due  to  inertia.  Therefore, 
the  trajectory  of  the  particle,  such  as  PQ',  deviates  from  the 
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streamline.  The  Inertia  varies  with  particle  diameter  and 
density.  Therefore,  the  extent  of  deviation  is  also  different. 

Some  particles  impact  the  plate  and  some  flow  through  the  hole  in 
the  next  stage. 

Assumption  (1):  The  flow  is  incompressible  and  inviscid. 
Moreover,  the  effect  of  gravity  is  neglected. 

Assumption  (2):  The  particle  velocity  has  already  caught  up 
with  the  flow  velocity  before  reaching  the  jet  outlet.  It  is  in 
equilibrium.  Particles  are  uniformly  distributed  in  the  flow. 

At  the  jet  outlet,  the  flow  velocity  is  uniform. 

Assumption  (3):  The  effect  of  diffusion  is  neglected. 
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Figure  1 


Assumption  (4):  The  particle  content  is  very  small  and  the 


presence  of  particles  does  not  affect  the  flow  field. 
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Parameters  affecting  the  motion  of  particles  are: 

D  -  aperture  diameter,  v^-  jet  velocity, 

p  -  gas  density,  u  -aerodynamic  viscosity  coefficient, 

Pp-  density  of  particle  material,  d  -  particle  diameter, 

2s-  collector  plate  diameter,  h  -  impact  distance, 

g  -  gravitational  acceleration,  a  -  speed  of  sound, 

A  -  mean  free  path  of  gas  molecule, 

xp0 -  initial  value  of  radial  coordinate  of  the  particle 
Xp(i.e.,  at  the  outlet  of  the  jet), 
t  -  time. 

Let  A  represent  a  function  of  the  velocity  components  of  the 
particle  u  and  v  or  the  coordinates  of  the  particle  x  and  y  . 

r  r  r  r 

Its  general  form  is 

A  -  h(D,  vm,  p,  ft,  Of,  d,  if,  t ,  ;  1 >  **•  0 

When  studying  the  trajectory  of  a  particle,  then 

.  * 

*!  "  f(.D,  *m,  P,P,  Pf,  d.  If,  H,  g,  •»  *#•»  y») 

The  parameters  in  the  above  formula  are  rendered  non-dimensional 
by  using  D,  vm  and  p.  By  taking  assumptions  (1)  and  (2)  into 
account,  we  get: 

X,-f( 9,  ®±,S,H,  K„  X„,  y,J 


where®  =  pDv^/y  (Reynolds  number) 


Kn=\/d  (Kenuzhen  Constant) 
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S  -  2 t/D, 

X,  -  2 x,/D, 
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H  -  IkID, 

-  2 xJD, 
Y,  -  lyP!D. 


In  the  limiting  case  that  both  H  and  S  are  very  large,  the  above 
formula  is  simplified  as: 

X,  -  G  (s„  Kn,  X *,  y ,) 

The  Impact  point  of  the  particle  on  the  collector  plate  is 

Xft  “  G  (^Stt X*  K»t  Xptt  o) 


When  Xpc-  •»  V  "<XP#)cr*  Then 

(Xja  -  Q  (s„  ® ±,  Ku'j 

Because  of  Assumption  (2),  the  impact  coefficient  P=(Xp0  )^.r*  It  /465 
represents  the  ratio  of  the  number  of  particles  at  a  certain 
diameter  collected  on  the  collector  plate  to  the  total  number  of 
such  particles.  In  sampler  studies,  it  is  called  the  relative 
collecting  probability^^.  Computation  shows  that  the  solution 
obtained  with  a  very  large  H  is  still  appropriate  at  H=1.5.  In 
the  following,  functions  G  and  Q  are  determined  numerically  based 
on  fluid  dynamic  equations. 


II.  Basic  Equations  and  Boundary  Conditions  of  the  Flow  Field 
In  an  axi-symmetric  coordinate  system,  the  velocity 
potential  f  of  an  inviscid  and  incompressible  flow  satisfies  the 


following  equations: 


&2L  +  +  1 

Ar1  dy*  *  3* 


0 


*  t  ♦ 

•."WV 


25 


where  x  is  the  radial  coordinate  and  y  is  the  axial  coordinate. 
The  subscript  g  indicates  gas.  CD  and  AB  are  equi-potential 
surfaces.  DO  and  OA  are  zero  flow  lines.  CB  is  the  free 
boundary  of  the  jet  (see  Figure  2).  The  boundary  conditions  are: 
9  *  constant  (i.e.,  Ug-0 ,  Vg=-vo8)  on  CD 
9  '«  constant  (i.e.,  Ug=v-f  Vg=0)  on  AB 
9  *  constant  (i.e.,  nx .  Vf=0),  x.Vi=vm  on  BC 
9*0  (i.e.,  nt  .  V9*0)  on  DO  and  OA. 
where  9  is  the  flow  function,  n!  is  the  unit  vector  in  the  normal 
direction  and  is  the  unit  vector  in  the  tangential  direction. 


Figure  2 


The  velocity  components  are: 
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(2.6) 


where  n  and  t  are  the  unit  vectors  in  the  normal  and  tangential 
direction  on  the  XY  plane. 

This  is  a  LaPlace  equation  with  unknown  boundary  BC.  By 
using  the  character is itcs  of  the  harmonic  function,  the  Green 
function  is  introduced  to  convert  the  above  differential  equation 
to  an  integral  equation.  Under  specific  conditions,  it  can  be 
integrated.  This  is  the  so-called  point  source  function  method 
in  fluid  dynamics. 

-  The  unknown  boundary  BC  can  be  assumed  to  be  a  known  curve. 
The  velocity  potential  (t)  of  various  points  on  BC  and  the 
potential  velocity  of  points  on  DC  and  AB  are  calculated. 
Furthermore,  the  velocity  potential  of  points  on  OA  is  also 
calculated.  Based  on  the  integral  expression  of  velocity 
potential,  the  velocity  potential  *a(t)  of  points  on  BC  can  be 
calculated.  The  shape  BC  can  be  repeatedly  adjusted  until  (t)*» 

In  our  calculation,  we  chose  Y(C)*Y(D)*3.  In  this  case,  at 


X(A)«x(B)»4,  the  flow  velocity  is  already  very  uniform. 
Moreover,  the  non-dimensional  velocity  *1.  The  result  of  this 
calculation  is  shown  in  Figure  3. 
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(3.3) 

(3.4) 


Wg  and  Wp  are  the  velocity  vector  of  the  primary  flow  and  the 
velocity  vector  of  the  particle,  respectively, 
rp  is  the  particle  position  vector, 
t  is  the  time. 

For  a  small  particle,  the  drag  coefficient  CD  can  be 
expressed  as 

(3.5) 


C,-7-Kh)/» 


where  Re*  d.p.w __/u  which  is  the  Reynolds  number  of  the  particle 

-  '  m  Or 

based  on  using  the  particle  diameter  d  as  the  characteristic 
length  and  relative  velocity  wgp  as  the  characteristic  velocity. 

w  -  1  +  2.46x/d  (3.6) 

where  w  is  a  correction  factor  for  the  dilute  gas  effect  and  x  is 
the  mean  free  path  of  the  gas  molecule.  Obviously,  when  X<<d, 
w*1 . 

In  this  calculation,  we  chose  f(Re)*1  and  f  (Re)=>1+1/6Re^^. 
Therefore , 
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and 
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(3.7a) 


(3.7b) 

Equations  (3.1)  and  (3.2)  are  made  non-dimensional  by  using 
D/2  and  vm  as  the  characteristic  quantities. 
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(3.11) 


w(g)  (Stokes  number) 
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R,  -  2  rf(D 


The  relation  between  R  and  Re  is: 


Re  «  R  d/D  |W, 


(3.17) 


(3.18) 


Ty  is  the  relaxation  time  used  to  judge  the  lagging  of  the 
particle  velocity  variation  behind  the  flow  velocity  variation. 
(S)and  St  can  be  considered  as  non-dimensional  relaxation  times. 
From  equation  (3.8)  we  can  see  that  when  St  is  very  small,  i.e. 
when  the  relaxation  time  is  very  short,  the  velocity  of  the 
particle  and  that  of  the  flow  are  in  equilibrium.  Otherwise,  it 
is  not  in  equilibrium. 

When  cartesian  coordinates  X  and  Y  are  used  on  the  azimuthal 


plane,  the  equations  become 


(3.19) 


s 


when  T-0,  Up«0,  Vp—  1,  Xp*Xp0  ,  Yp-Y(C)=H  (3.23) 

U  .  and  U.,  V„  are  the  X,  Y  components  of  W„  and  W  . 

P  P  8  8  P  g 

respectively.  Xp  and  Yp  are  the  components  of  Rp. 

On  the  meridian  plane,  if  the  potential  function  *  and  flow 
function  T  are  expressed  in  an  orthogonal  coordinate  system,  then 


the  equations  become 

dW„ 

' 7T 
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When  T*0,  W  *1,  W  *0 ,  •  «*(C)=*A,  f  =f  . 

pt  pn  p  p  p* 


(3.27) 


(3.28) 


where  ,  W  ,  ♦  and  (-Y  )  are  the  components  of  W  and  R  in 
Pt  pn’  p  p  r  p  p 

the  direction  of  t  and  n,  respectively,  t  is  the  unit  vector  in 
the  tangential  direction  of  the  streamline  pointing  in  the 
direction  of  increasing  ♦  .  n  is  the  unit  vector  in  the 
streamline  direction  pointing  in  the  direction  of  decreasing  f. 
a  is  the  angle  of  rotation  from  the  x-axis  counterclockwise  to  t- 
direction.  The  problem  is  that  a<0. 


2.  Integration  of  Equation  of  Motion  of  Particles 
(1)  Method  of  Integration 


For  a  very  small  particle  (e.g.,  St<0.05),  because  its 

trajectory  of  motion  almost  coincides  with  the  streamline,  it  is 

2 

possible  to  employ  a  perturbation  method.  When  (2St)  is  very 
small  (in  general  f(Re)**1)  and  d  Up/dT* is  not  very  large,  we  get 
the  following  from  equation  (3.19) 
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A  similar  expression  can  be  obtained  for  Vp. 

-  If  aT  (step  length)»0.15,  when  0.01<St<0.075 ,  equations 

(3.19)  -  (3.22)  can  be  integrated  using  the  Treanor  method^"®. 

In  our  computation,  it  was  proven  that  the  results  obtained  with 

Treanor' s  method  are  in  total  agreement  with  those  obtained  with 

the  perturbation  method  at  St=0.01  and  aT»0.15. 

When  St>0.075-,  the  R-K  method  can  be  used.' 

For  very  small  particles,  a  larger  step  can  be  used  in  the 

perturbation  method.  In  this  case,  it  is  more  appropriate  to  use 

the  •-V  coordinate  because  is  very  small  and  W  does  not 

pn  J  pr 

vary  significantly. 

Thus,  a  series  of  particle  trajectories  can  be  calculated 
corresponding  to  any  given  St  and  R  d/D.  Taking  the  requirements 
of  sampler  design  into  account,  we  used  n*  to  replace  R  d/D: 

n*  =  (R  d/D)/St  (A) 


(2)  Selection  of  Step  Length 

In  addition  to  considering  requirements  of  satisfying  the  / 469 
stability  of  the  difference  scheme  and  the  accuracy  of  the 
computation,  we  should  try  to  shorten  the  computing  time  to  the 
extent  possible  in  order  to  make  it  practical.  In  this 
calculation,  we  chose  aT=0.15.  The  average  time  to  calculate  a 
trajectory  is  60  seconds  (on  a  Model  FILEX-512  general  purpose 
computer).  Table  1  shows  the  Y  value  at  x=4  on  each  trajectory 
when  S^=0  (in  this  case,  trajectories  are  streamlines).  These 
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values  were  compared  to  those  corresponding  to  Y=1/8Xp,  and 
were  found  in  good  agreement.  This  demonstrated  that  the 
required  accuracy  is  satisfied  in  the  calculation. 


(3)  Several  Measures  to  Cut  Down  Computing  Time 

It  took  more  than  10  minutes  to  calculate  each  trajectory 
using  the  initial  program.  The  following  measures  were  taken  to 
reduce  it  to  less  than  1  minute. 

(a)  Frequently  used  elliptical  Integrals  are  calculated  by 
series  expansion  to  drastically  save  time  and  obtain  good 
accuracy. 

(b)  Computation  of  Flow  Velocity  and  its  Partial 
Derivatives 

For  Y  >0.03,  the  flow  velocity  is  calculated  based  on  its 
integral  expression.  For  Y^0.03,  the  flow  velocity  at  any  point 
is  calculated  by  series  expansion.  Not  only  time  is  saved  but 
also  accuracy  is  improved.  Our  computation  showed  that  the 
series  expansion  method  can  be  extended  to  Y<0.3. 

As  for  the  calculation  of  various  coefficients  in  the  series 
expansion  equation,  i.e.,  U  ,  aU  /dX,  a*U  /aX*  ,  a*  U  /aX*  at 

o  o  o  o 

various  points  on -the  x-axis,  the  first  two  are- calculated  by 
using  the  integral  expressions  and  the  latter  two  are  calculated 
first  by  using  a  sample  value  fitting  method  with  discrete  values 

of  U  and  aU  /aX  and  then  by  differentiation. 

8  g 

We  must  point  out  here  that  even  though  a  particle  he 
already  reached  Y=0.03,  we  still  cannot  get  the  approximate  value 
at  Y=0  by  extrapolation.  When  Y<  0.03,  the  particle  trajectory 
bends  significantly.  Furthermore,  for  particles  in  the  range  of 
St=0. 1—0.2  (which  is  the  range  of  our  concern),  is  very  small 
when  Y»0.03.  Therefore,  extrapolation  is  not  reliable. 


v.iVvvV. 
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(c)  The  computation  started  from  Y=1.5  due  to  the  zt  that 


the  flow  field  is  very  uniform  at  Y»2  to  Y=3,  the  plate  has  no 
effect  on  the  jet,  and  the  particle  trajectory  is  straight. 

IV.  Analysis  of  Calculated  Results 
Figure  3  shows  the  streamlines.  Figure  4  shows  the 
calculated  particle  trajectories  for  St=0.12.  If  there  is  an 
arrow  on  a  trajectory  near  the  boundary  OA,  then  this  trajectory 
does  not  intersect  with  OA  within  X<3.5.  Figure  5  shows  the 
effect  of  0*  on  the  particle  trajectory.  Our  computation  shows 
that 

/470 


Figure  4.  Particle  Trajectories 


1.  jet  boundary 


(1)  Within  the  range  of  X<3.5,  almost  all  the  particles 


with  St<0.11  cannot  reach  the  boundary  OA  unless  Xp„  is  very 
small.  Almost  all  particles  with  St>0.17  can  reach  the  OA 
boundary.  Whether  particles  in  the  range  of  0.11 <S^0 .17  can 

reach  the  boundary  OA  at  X<3.5  depends  on  the  values  of  Xp0  and  0 
* 


(2)  Impact  points  are  basically  within  X<2.0.  At  X>2.0, 

impact  points  are  obviously  frequent.  This  effect  can  be 

explained  from  equation  (3.24)  and  (3.25). 

The  impact  points  of  particles  with  large  St  are  within 

X=2.0.  Let  us  now  discuss  particles  with  St<0.2.  In  most  areas 

of  the  flow  field,  | aa/3T | << | 3a/3* | .  For  smaller  particles,  Wpn<< 

W^  .  Thus 
P* 


doLmm}_w  da  W„  da  W,r  da 

jt  p,  "  &p  p,  dr  p,  d<P 


W,r/R' 


(4.1) 


where  is  the  radius  of  curvature  of  the  streamline.  By 
substituting  equation  (4.1)  into  (3.25),  we  get 

dW  M  t»i  Ip  j(.Re )  nr 

-£r-W„lRw  2S'  w„  (4.2) 

By  omitting  the  second  term  on  the  right  of  equation  (3.24),  we 

8  -  liEsl  (yt  -  w„) 

n  251  (4.3) 


The  first  term  on  the  right  side  of  equation  (4.2)  represents  the 
inertia  centrifugal  force  which  is  the  driving  force  to  increase 
Wpn.  The  second  term  on  the  right  represents  the  aerodynamic 
drag  which  is  the  damping  force  to  decrease  Wpn>  When  the 
centrifugal  force  can  be  neglected,  Wpn  decays  exponentially. 
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The  decay  constant  is  S^Cf(Re)al).  The  smaller  the  particle  size 
is,  the  faster  the  decay  becomes.  For  a  microscopic  gas  group,  S 
=0.  Thus,  VJ  =0.  For  a  particle  (S^>0),  W„  first  gradually 
Increases  from  zero  and  then  gradually  decays  to  zero.  is 

always  decreasing.  When  *p=0,  the  particle  reaches  the  collector 
plate. 


Figure  5.  Particle  Trajectories 


Siigl 

»  ! 

When  T-T  ,  W  >>2S.W*  / r 

a*  pn  t  p-r  ¥ 

(4.4) 

VO 

From  (4.2)  we  get 

v  ■%  .*1 
* 

w fm  —  (H',.)*«p  [— 

KRO  (t  —  t.)1 

2 s,  -1 

(4.5) 

Substituting  it  into  equation 

(3.27) ,  we  get 

[■  -  ■* [-  (r  TJ] 

}  (4.6) 

,■  V/, 

The  subscript  a  represents  the 

value  at  T=T .  When  (T  )  ,  (W  ) 
a  pa  pn  a 

.\\K* 

and  St  are  given,  we  can  find 

the  value  of  (T-T  ) 

a 

corresponding 

to  (see  Table  2).  When  equation  (4.4)  is  satisfied,  usually 

W^«U  «1.  Therefore,  the  second  row  in  Table  2  corresponds  to  (X _ 

pt  p  ’  r  pc 

-Xpa)/2St.  It  is  obvious  from  Table  2  that  the  impact  points  of 
a  certain  particle  size  are  spaced  less  densely  away  from  the 
center  on  the  collector  plate.  This  point  was  also 
experimentally  verified. 


4-  •* 


Table  2 

C  * 


o.i 

0.7 

0.1 

0.9 

0.92S 

0.95 

0.97 

0.99 

1 

(T  - 

0.91 

1.20 

l.tl 

2.3 

2.39 

3.00 

3.30 

4.60 

GO 

Figure  6.  The  P(St,  n*)  Curve 

1.  not  corrected  for  viscosity 

2.  corrected  for  viscosity 

The  above  discussion  indicates  that  the  results  of  an 
infinitely  large  plate  (i.e.,  S  is  very  large)  dan  be  applied  to 
a  finite  plate  as  long  as  the  edge  effect  on  the  X<2.0  region  can 
be  neglected. 

(3)  Corresponding  to  a  specific  St  and  0*,  there  is  a 

certain  X  «  (which  is  de'noted  as  (X  »  )  )  which  corresponds  to  X^ 

p  p°  cr  r  p 

-3.5. 

From  p  =  (Xp0  )^r  (4.7) 

we  can  find  the  impact  coefficient  p(St,0*)  corresponding  to  that 
particular  S  and  0*. 


V.  Effect  of  Viscosity 

When  an  impact  type  collector  is  used  to  collect  particles 
larger  than  1>u  in  diameter,  the  Ma  number  is  usually  less  than 
0.2  and  Dg/v«  is  less  than  0.3.  The  flow  field  is  a  low  Reynolds 
number  laminar  flow  in  which  the  effect  of  diffusion  can  be 
neglected.  The  accumulation  of  a  few  particles  in  the  atmosphere 
will  only  affect  the  flow  in  the  boundary  layer.  It  has  little 
effect  on  the  inviscid  flow  field.  Therefore,  with  the  exception 
of  neglecting  the  effect  of  viscosity,  other  assumptions  made  in 
performing  the  calculation  are  reasonable. 

1.  Analysis  of  Viscosity  Effect 

-The  effect  of  viscosity  is  primarily  exhibited  in  three 

[4] 

areas  : 

(1)  Due  to  viscosity,  the  jet  outlet  velocity  is  non- 
uniform. 

(2)  There  is  momentum  exchange  across  the  free  boundary  BC. 
Momentum  is  transferred  from  the  inside  to  the  outside  of  BC. 

(3)  There  is  a  boundary  layer  at  the  wall  OA  of  the  plate 
which  makes  Ug  near  OA  smaller  than  the  value  in  an  inviscid 
flow.  Due  to  the  presence  of  the  solid  wall,  the  drag  on  the 
particle  is  also  affected. 

With  regard  to  the  effect  of  the  first  point,  the  following 
correction  is  made.  Let  us  assume  that  the  distribution  of 
velocity  at  the  outlet  of  the  jet  is: 


!•  the  velocity  at  the  center  of  the  jet  outlet,  v  is 
•  max  00 

the  mean  velocity  at  the  jet  outlet.  Hence 
St/St=1.264[1-(Xp0)cr]1/6 

where  is  defined  as  t.  /D  and  S„  is  defined  as  t„.v _ ,/D. 

t  V00  l  v  ®raax 

If  we  assume  that  the  particle  trajectory  passing  through  a  point 
at  the  outlet  only  depends  on  the  local  St  value  and  is  not 
related  to  the  fact  whether  the  jet  is  uniform  or  not,  then  the 
impact  coefficient  P  is 

P  -  f.  '  2xXdXl*Sm 

-  1  -  l4  11  -(XM”  ”  T  11  " 

o  o 
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With  regard  to  the  effect  of  item  (2),  because  the  flow 
velocity  distribution  near  the  boundary  (jet  mixing  area)  is 
affected  by  the  momentum  exchange  across  the  free  boundary,  the 
particle  trajectory  is  also  altered.  The  half-width  b  of  the 
mixing  area  at  H=1.5D  was  estimated  by  using  the  semi-inif inite 

[A] 

free  jet  formula  (choosing  velocity  ratio  =0.95  as  the  jet 
boundary).  It  was  found,  that  b/D/2  is  in  the  range  of  0.22-0.093 
which  shows  that  the  jet  boundary  only  affects  the  P>0.65-0.75 
portion  of  the  P(§t,fl*)  curve. 

With  regard  to  the  effect  of  item  (3),  due  to  the  presence 
of  the  boundary  layer  near  the  wall,  the  transverse  velocity  near 
the  wall  is  lower  than  that  of  an  inviscid  flow,  which  affects 
the  flow  of  particles  near  the  wall.  For  a  collector  operating 


at  a  flow  rate  of  2  liter /min,  the  Reynolds  number  Re  ranges 

from  626  (first  stage)  to  2606  (sixth  stage).  It  is  estimated 

that  the  ratio  of  the  boundary  layer  thickness  60  to  the  radius 

of  the  jet  D/2  ranges  from  0.176  (first  stage)  to  0.080  (sixth 

stage).  It  is  very  small  compared  to  H=3.  In  this  thin  layer, 

the  streamline  is  already  very  straight.  After  a  particle  enters 

this  thin  layer,  its  normal  direction  velocity  decays  very 

rapidly.  Tp  will  not  drop  significantly  any  further.  Therefore, 

the  effect  of  the  boundary  layer  is  mainly  causing  the  tangential 

velocity  of  the  particle  to  drop  which  consequently  leads  to  an 

increase  of  time  of  particle  motion.  It  has  little  effect  on  the 

1 ) 

shape  of  the  particle  trajectory  . 

2.  Corrected  Results 

The  viscosity  corrected  impact  coefficient  P(§t,Q*)  curve  is 
shown  as  the  solid  curve  in  Figure  6.  When  n*  varies  from  0  to 
250,  St,0  (the  value  of  St  corresponding  to  P=50%)  varies  from 
0.122  to  0.137.  Corresponding  to  our  calculating  conditions, 
Marple  and  Wileke^  got  St,  0  =0 . 1 165-0. 1 12^  when  O*=0.  R. 
Wiedemann^  choseN/  St,  0  =0.3433^  and  got  St,  0  =0. 1 1765.  In  some 

foreign  sampler  designs,  S^,,  =0.144  was  chosen.  These  are  in 
agreement  with  our  computation. 

The  calculated  P(§t,fl*)  curve  was  used  to  design  a  multi¬ 
stage  impact  sampler.  After  evaluation,  it  was  found  that  the 
performance  of  the  sampler  agreed  well  with  theoretical 
calculation.  Design  specifications  are  met^2^. 


1)  More  profound  investigation  not  only  requires  the 
consideration  of  the  inhomogeneity  of  velocity  in  the  boundary 
layer  but  also  must  include,  in  addition  to  the  drag  in  the 
direction  of  motion,  the  transversal  force  exerted  on  the 
particle  by  the  flow. 

2)  In  reference  [5],  Stjc,,/2  corresponds  to  St,0  in  this  work. 

3)  In  reference  [6],  Stk  corresponds  to  St, 0  in  this  work. 

4)  In  reference  C  7] ,  <A.  corresponds  to  St,0  in  this  work. 

VI.  Conclusions 

Uhen  there  are  fewer  particles  in  the  jet,  they  do  not 
affect  the  flow  field  significantly.  Therefore,  the  computation 
of  the  particle  trajectory  and  the  flow  field  can  be  separately 
discussed.  Under  the  inviscid  and  incompressible  conditions,  the  / 474 
flow  field  is  calculated  by  solving  an  unsteady  boundary  Laplace 
equation.  A  "point  source  function"  method  is  used  to  convert 
the  differential  equation  to  integrals.  The  shape  of  the  unknown 
boundary  is  repeatedly  adjusted  to  reach  the  correct  position. 

The  particle  trajectory  is  calculated  by  solving  a  series  of 
regular  differential  equations.  Based  on  the  magnitude  of  St, 
either  the  perturbation  method,  Treanor  method  or  R-K  method  can 
be  used.  Our  computation  shows  that  the  impact  points  of  the 
particles  on  the  plate  are  primarily  concentrated  in  the  range  X< 

2.  The  calculation  made  for  very  large  H  is  still  valid  at 
H*1.5.  After  taking  the  effect  of  viscosity  which  causes  non¬ 
uniformity  of  the  velocity  at  the  jet  outlet  into  account,  a 


velocity  distribution  va>/veB=1.264(1-2x/D)  '  is  used  to  correct 
the  results.  A  P(§t,  a*)  curve  is  obtained.  A  sampler  was 
designed  according  to  that  curve  and  its  performance  met  the 

design  specifications^. 

This  work  was  completed  under  the  guidance  and  assistance  of 
Professor  Wu  Chenkang.  Professor  Ban  Yinggui  also  provided 
guidance  and  assistance.  Comrade  Li  Jiachun  provided  many 
valuable  opinions  on  the  manuscript.  The  author  wishes  to  thank 
them  for  their  efforts. 
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Abstract 

The  calculation  of  the  flow  field  of  a  uniform  circular  jet  with  particles  impacting 
upon  an  infinite  plate  perpendicularly  is  presented.  The  fluid  is  assumed  to  be  inviscid 
and  incompressible.  The  trajectories  of  small  spherical  particles  carried  in  the  jet  are 

also  calculated  by  choosing  either  ~  0r  —  (l+Bew)  as  drag  coefficient. 

Be  Re 

Assuming  that  the  particles  are  uniformly  distributed  initially  in  the  jet.  we  obtain 
the  curves  of  impact  efficiency  (known  aa  collection  probability  in  the  study  of  sam¬ 
plers).  The  basis  upon  which  the  assumptions  are  made  in  the  calculation  is  discussed. 
Viscous  effects  are  analyzed,  and  then  some  corrections  of  the  curves  P(S„  £1*)  consi¬ 
dering  these  effects  are  made.  The  results  of  calculation  show  that  the  unpact  points 
of  the  particles  on  the  plate  mainly  concentrate  in  the  zone  of  x<2.  The  results  of 
calculations  baaed  on  very  large  S  is  valid  for  H  >1.5.  An  appraisal  for  the  sampler 
designed  by  using  P(S„  fi#)  curvet  shows  that  expected  performance  has  been  realized. 
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A  Non-local  Elastic  Plastic  Continuum  Model  and  the  Distress  ' 
Distribution  Near  a  Cracked  Tip 
Yu  Jilin  (University  of  Science  and  Technology  of  China) 
and  Zheng  Zhemin  (Institute  of  Mechanics,  Academia  Sinica) 

Abstract 

A  non-local  elastic  plastic  continuum  model  is  presented. 

In  this  model,  the  relation  between  stress  and  elastic  strain  is 
non-linear  and  plastic  strain  is  related  to  the  history  of  total 
strain.  With  regard  to  the  deformation  theory,  it  is  assumed 
that  the  plastic  strain  tensor  is  proportional  to  the  total 
strain  deviation  tensor.  The  proportionality  factor  is  a  scalar 
function  of  the  total  effective  strain.  This  model  was  used  to 
analyze  the  stress  field  at  the  tip  of  a  power-law  hardening 
material  with  a  tensile  crack.  Based  on  the  results  of  .RR 
asymptotic  solution  of  the  tensile  cracked  tip  obtained  in 
classical  fluid  dynamics,  expressions  for  the  distribution  of 
tensile  stress  in  front  of  the  crack  and  the  maximum  tensile 
stress  were  derived  through  one-dimensional  simplification.  They 
showed  that  the  J  integral  criterion  might  be  obtained  from  the 
maximum  tensile  stress  criterion.  Existing  experimental  data  were 
used  to  calculate  the  maximum  stress  near  the  tip  of  several 
steel  materials  as  the  crack  begins  to  propagate.  It  was 
discovered  that  its  order  of  magnitude  is  close  to  that  of  the 
cohesive  strength  of  the  lattice.  The  results  obtained  are 
useful  for  the  understanding  of  the  physical  mechanism  of  the 
fracture  process  of  the  material. 


I.  Introduction 


In  the  recent  two  decades,  based  on  classical  elastic 
mechanics  and  elastic  plastic  mechanics,  fracture  mechanics  has 
been  developed  at  a  rapid  pace.  It  serves  as  a  new  theoretical 
basis  for  the  safety  design  of  engineering  components,  estimation 
of  useful  life  and  evaluation  of  the  performance  of  engineering 
materials.  It  is  very  successful  in  engineering  applications. 

Classical  fracture  mechanics  was  developed  based  on  the 
equilibrium  criterion  which  was  proposed  by  Griffith  and  extended 
by  Orowan.  The  concept  of  J  integral  is  also  based  on  energy 
analysis.  This  type  of  energy  criterion  avoids  the  physical 
mechanism. of  the  fracture  process,  i.e.,  the  stress  and  strain 
conditions  near  the  cracked  tip.  In  reality,  according  to 
classical  continuum  mechanics,  there  is  a  stress  singularity  at 
the  tip.  The  widely  used  concept  of  stress  strength  factor  and 
the  HRR  analysis  of  power-law  hardening  materials  recognize  this 
singularity.  However,  if  this  singularity  exists,  the  cracked 
body  cannot  sustain  any  load.  This  contradiction  is  one  of  the 
major  deficiencies  of  classical  fracture  mechanics. 

Many  attempts  were  made  to  eliminate  such  stress 

rii 

singularity,  including  the  linear  yield  band  model  of  Dugdal  , 

T  21 

radius  of  curvature  correction  made  by  Neuber  ,  dislocation 

C31 

model  of  Bilby  et  al  and  super  dislocation  model  of  Atkinson 
et  al^^.  These  models  mostly  involve  local  correction  on  the 
basis  of  classical  theory.  They  are  useful  for  certain  practical 
problems.  However,  the  physical  basis  is  ambiguous. 


Some  people1-  *  J 


believe  that  there  Is  always  a  radius  of 


curvature  at  the  tip  crack.  The  ideal  sharp  crack  does  not  exist. 
Microscopically  (on  the  order  of  um) ,  it  may  be  the  case. 

However,  the  fracture  process  is  essentially  the  destruction  of 
the  atomic  bond.  On  a  finer  scale,  plasticity  does  not  flow 
uniformly.  There  are  dislocation  cells  when  the  strain  level  is 
high.  There  is  a  theory*-^  which  proves  that  in  most  metals,  with 
the  exception  of  face-centered  metals  and  alkali  metals,  the 
sharp  crack  on  the  atomic  scale  will  not  become  dull.  Iron  is  in 
the  middle.  Therefore,  despite  the  possible  plastic  deformation 
near  the  tip,  the  crack 

Manuscript  received  on  January  16,  1984. 

may  still  remain  sharp  when  the  accuracy  is  on  the  atomic  scale. /486 
In  reality,  a  material  is  composed  of  discrete  atoms.  It 
has  a  complex  internal  structure.  Physical  quantities 
corresponding  to  the  continuum  field,  such  as  displacement, 
strain,  stress,  etc.,  can  only  be  established  on  the  basis  of 
local  averages.  When  the  characteristic  scale  of  the  physical 
phenomenon  under  consideration  is  comparable  to  that  of  the 
internal  structure  of  the  material,  classical  continuum  mechanics 
will  encounter  difficulties.  A  more  rational  physical  structure 
theory  must  take  the  internal  structure  of  the  material  into 
account.  To  this  end,  since  the  60's,  various  continuum  theories 
which  take  the  micro-structure  of  the  material  in  consideration 
were  developed.  In  treating  problems  relating  to  cracked  tip  and 
its  stress  singularity,  the  non-local  theory  developed  by  Eringen 
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et  al  *  has  been  successfully  applied. 

The  non-local  theory  considers  that  the  interaction  between 
atoms  is  a  long  range  force.  Therefore,  the  stress  on  a  point  is 
not  only  related  to  the  strain  and  strain  history  at  that  point, 
but  also  related  to  the  strain  and  strain  history  at  other  points 
in  the  object.  In  other  words,  the  stress  at  a  point  is  a 
general  function  of  the  strain  field  of  the  entire  object  and  its 
history. 

Eringen  and  his  colleagues  used  the  non-local  theory  to 

study  the  stress  field  at  the  tip  of  a  brittle  single  crystal 

material.  Their  results  showed  that  the  stress  singularity  at 

the  cracked  tip  did  not  exist.  The  maximum  stress  appeared  at  a 

small  distance  in  front  of  the  cracked  tip.  In  addition,  the 

theoretical  cohesive  stress  determined  is  in  agreement  with  those 

obtained  by  using  atomic  theory  and  by  experimental  prediction^10^  •  ' 

means  that  the  cracked  body  also  meets  the  maximum  stress 

fracture  criterion.  Based  on  the  non-local  elastic  theory,  it  is 

not  necessary  to  introduce  the  surface  energy  which  has  no  clear 

physical  significance  in  classical  theory  to  directly  derive  the 

Griffth  criterion  from  the  maximum  stress  criterion.  Therefore,  the 
physical  significance  of  the  finding  is  very  high. 

In  this  work,  the  non-local  theory  is  used  to  study  the 
cracking  problem  in  elastic  plastic  materials.  When  plastic 
deformation  is  invlolved,  the  state  of  stress  is  related  to  the 
strain  history.  Let  us  consider  a  metallic  lattice  (single 
crystalline  or  polycrystalline).  Under  an  external  load,  if 
there  is  only  elastic  deformation,  the  spacing  between 
neighboring  atoms  in  the  material  deviates  only  very  slightly 
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relative  to  its  equilibrium  value.  Once  the  external  load  is 
removed,  these  atoms  can  still  return  to  their  original  states. 
However,  plastic  flow  is  related  to  the  motion  and  increase  of 
dislocations  in  the  crystalline  material.  In  other  words,  the 
plastic  deformation  of  material  signifies  permanent  changes  of 
the  atomic  arrangement,  i.e.,  relative  positions  of  atoms,  in  the 
material.  As  a  plastic  flow  develops,  there  is  a  new  equilibrium 
in  every  instance.  Due  to  the  presence  of  the  external  load, 
there  is  a  small  deviation  of  the  atomic  spacing  away  from  the 
new  equilibrium  state.  Based  on  this  physical  picture,  stress  is 
corresponding  to  the  deviation  of  atomic  spacing  away  from  the 
equilibrium  state  at  that  instance,  i.e.,  the  elastic  portion  of 
strain.  Considering  the  fact  that  the  interaction  between  atoms 
is  a  long  range  force,  stress  and  elastic  strain  should  be 
described  by  a  non-local  relationship.  However,  plastic  strain 
corresponds  to  the  permanent  change  of  equilibrium  atomic 
arrangement.  It  should  be  related  to  the  total  'Strain  history. 
Based  on  the  small  deformation  theory  and  plastic  deformation 
theory,  the  structure  and  basic  equations  for  this  type  of  power- 
law  hardening  non-local  elastic  plastic  material  were  established 
in  this  work.  Furthermore,  it  was  used  to  study  the  stress  field 
near  the  tip  with  a  tensile  crack.  The  HRR  singularity  solution 
of  a  tensile  crack  near  the  tip  obtained  by  classical  fracture 
mechanics  was  used  to  derive  the  expression  for  the  maximum 
stress  near  the  tip  in  small  scale  yield  under  one-dimensional 
simplification.  It  has  been  proven  that  the  critical  J  integral 
criterion  in  classical  fracture  mechanics  can  be  obtained  from 


the  maximum  stress  criterion  on  the  non-local  elastic  plastic 
theory.  Existing  experimental  datawere  used  to  calculate  the 
critical  maximum  stress  values  near  the  cracked  tips  of  several 
steel  materials.  It  was  discovered  that  it  is  of  the  same  order 
of  magnitude  as  the  internal  cohesive  stress  of  the  lattice. 
Finally,  the  physical  meanings  of  this  theory  and  its  results  are 
discussed. 


II.  Basic  Equations 

In  this  work,  the  plastic  structure  is  described  in  the 
strain  space  because,  as  described  in  the  previous  section,  that 
plastic  deformation  is  related  to  the  total  strain,  rather  than 
stress . 

According  to  the  small  strain  theory,  strain  e  is  related  to 
displacement  u  linearly,  i.e. 

«tj  -  1/2  <Uljj.  Ujl)  (1) 
where  the  subscript  following  the  comma  represents  the  partial 
derivative  with  respect  to  the  corresponding  coordinate. 

The  strain  is  divided  into  two  parts,  i.e.  elastic  strain 
and  plastic  strain: 


e  =  ee  + 
eij  cij  +  eij 


(2) 


where  the  superscripts  e  and  p  represent  elastic  and  plastic, 
respectively.  Let  us  assume  that  elastic  strain  ee  and  stress  t 
satisfy  a  non-local  linear  relation  [8] 

fcij  <*>-#*•  (  lx’  "xl  )ekk^x'  C  |x’  -x  |  )e®  j  (x' )  ]dV(x ' )  (3) 

where  is  Kronecher's  6,  and  a'  and  m'  are  non-local  modulus. 
All  repeating  subscripts  indicate  the  summation  with  respect  to 


all  Indices.  Let  us  assume  that  the  effect  of  factors  such  as 
the  Internal  characteristic  scale  change  of  the  materials  and 
microscopic  inhomogeneity  due  to  plastic  deformation  on  non-local 
modulus  can  be  neglected.  Then,  X'  and  y '  is  only  a  function  of 
position  |x'-x| . 

When  there  is  no  volumetric  force,  the  stress  equilibrium 
equation  is 

Hj.1  -°  (4> 

In  order  to  discuss  the  plastic  structure,  we  introduced 
strain  deviation  e  and  elastic  strain  deviation  ee 

i  i  (5) 

€H  -  •;/  —  —  ««*.!.  «<»  “  "  7  *u5ii 


as  well  as  effective  strain  t  and  effective  elastic  strain  e 


-e 


(6) 


Corresponding  to  the  classical  plastic  deformation  (full)  theory, 
in  the  strain  space  we  assume  that: 

1.  The  volumetric  strain  is  elastic,  i.e.,  e^=0.  Plastic 
strain  is  only  related  to  the  deviation  of  total  strain. 

2.  The  effective  elastic  strain  is  a  well-defined  function 
of  the  effective  total  strain. 

3.  The  plastic  strain  tensor  is  proportional  to  the  total 


strain  deviation  tensor.  The  proportionality  constant  is  the 
scalar  function  of  the  total  strain,  i.e. 
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Based  on  assumption  (2)  above,  this  relation  is  still  valid  for  a 
non-uniform  deformation  (in  this  case  °ij^*  From  equations 

(7)  and  (12)  we  can  solve  that 

{0  when  e<e„ 

n-1  (13) 

1-(e/c  )n  when  e>ey 

Thus,  equations  (1)-(7)  and  (13)  form  the  basic  equations  in 
the  non-local  elastic-plastic  deformation  theory  of  power-law 
hardening  materials  which  are  incompressible  and  are  deformed 
slightly.  It  is  only  applicable  to  simple  loading  and  near 
simple  loading  situations. 

III.  Stress  Field  Near  the  Tip  of  a  Tensile  Crack  Under  Plane 
Strain 

Let  us  assume  that  there  is  a  crack,  2a  in  length,  on  an 
infinite  plate.  The  surface  is  free.  At  infinity,  there  is  a 
uniform  tensile  stress  t<=°  perpendicular  to  the  crack  plane  (see 
Figure  1).  Near  the  cracked  tip,  plastic  strain  is  much  larger 
than  elastic  strain.  Therefore,  the  incompressibility  assumption 
is  approximately  valid.  If  we  are  limited  to  the  study  of  the 
Initial  propagation  of  a  crack,  deformation  theory  is  also 
applicable  to  power-law  hardening  materials.  When  the  dimension 
of  the  plastic  region  is  smaller  than  the  thickness  of  the  plate, 
it  can  be  approximately  considered  as  in  a  plane  strain  state. 
Because  of  symmetry,  we  are  only  required  to  study  the  upper 
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(14) 


The  problem  Is  reduced  to  solving  the  equations  in  two-dimensions 
under  the  following  boundary  conditions: 


*11**12"®*  t22=t»’  when  (xj+xj)  -* 

*12**22“®’  when  Xa*0»  |*k  |<  e 

8ui  /  dx»  =0 ,  ua  =0 ,  when  xa  =0 ,  |  x&  |  ^a 


The  last  condition  is  obtained  from  symmetry. 


Figure  1.  The  Tensile  Crack  Problem 
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Eringen  proved  that  the  non-local  modulus  can  be  /489 


expressed  as 


i' —  ia(|*’ —  *!)»  fi  — 


(16) 


where  a(|x'-x|)  is  a  non-negative  function  of  x1  which  has  the 
following  properties: 

1.  With  increasing  |x'-x|,  «(|x'-x|)  approaches  zero  very 
rapidly. 

2.  In  the  extremum  case  in  classical  elastic  theory, 
a(|x'-x|)  becomes  a  Dirac  6  function. 

3.  /ya( |x*-x| )dV(x* )=1 . 

By  comparing  to  the  ideal  lattice  model,  we  get 

,  t  (17) 

o  l*'-*l>* 

where  b  is  the  lattice  parameter  and  K=3/flba  (two-dimensional)  or 

K»1/b  (one-dimensional). 

Based  on  equation  (16),  the  equilibrium  equation  (4)  can  be 
rewritten  as  , 

.(I*-*|).*i(*W*')-<> 

]r  (18) 
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It  can  be  proven  that  if  the  effect  of  surface  tension  can  be 
neglected  the  necessary  and  sufficient  condition  for  the  above 
equation  to  be  satisfied  is: 

°ij , j  *°  <19) 

This  shows  that  the  control  equation  for  the  nominal  non-local 

stress  is  the  same  as  that  in  classical  theory.  Thus,  if  the 

effect  of  surface  tension  is  neglected,  for  a  given  displacement 

and  boundary  condition,  the  displacement  field  and  strain  field 

obtained  using  the  model  described  in  this  paper  are  identical  to 

those  obtained  based  on  the  classical  elastic  plastic  theory. 

However,  the  stress  field  must  be  calculated  using  equation  (14). 


The  above  conclusions  is  generally  not  valid  for  problems 


with  given  stress  boundary  conditions  or  mixed  boundary 

conditions.  If  classical  strain  field  is  used  to  calculate  the 

non-local  stress  field,  stress  boundary  conditions  can  only  be 

approximately  satisfied.  Despite  so,  because  the  microscopic 

scale  reflecting  the  non-local  effect  is  very  small  (atomic 

spacing),  actual  errors  only  occur  at  places  where  classical 

theoretical  stress  has  breakdowns  or  irregularities  (of  course, 

those  are  the  areas  of  concern) .  The  calculations  performed  by 
roi 

Eringen  on  a  non-local  elastic  crack  showed  that  this  error 
decreased  with  increasing  2a/b.  When  2a/b=40  (equivalent  to  a 
crack  length  around  O.lum),  the  maximum  stress  error  is  about 
10%.  The  actual  crack  is  much  larger.  The  error  will  be  less. 

For  power-law  hardening  materials,  based  on  classical  theory 
of  elastic  plastic  deformation,  the  stress  and  strain  at  the  tip 
have  r“n/1+n  and  r-1/1+n  singUiarities ,  respectively.  According 
to  reference  [12],- in  small  range  yield  cases,  the  major  terms  of 
stress  and  strain  fields  near  the  tip  can  be  expressed  as: 


where  v  Is  the  elastic  Polssonlt  of  the  material  and  In  o  is  a  constant 
related  to  the  hardening  index  n  (see  Figure  2).  cf  (e)  , . . . , 

(e)are  angular  distribution  functions  corresponding  to  °rr»***er0 
which  are  also  related  to  n.  The  angular  distribution  functions 
of  three  stress  values  at  different  n  are  shown  in  Figure  3. 


Figure  2.  I  vs.  n 


Figure  3.-  Stress  vs.  6  at  Various  Hardening  Indices  (from 
reference  [12]) 


If  equation  (21)  can  also  be  considered  as  an  acceptable 
approximate  solution  of  the  strain  field  in  non-local  elastic 
plastic  theory,  based  on  equation  (20)  we  can  find  the 
approximate  solution  of  the  non-local  stress  field  from  the 
following  equation: 


'a  “  \,  •(!*'  “  *1  >, ■<(*')*(*') 


(22) 


This  is  a  very  difficult  task  which  will  not  be  investigated 
further  in  this  paper.  _0ne  can  see  that  tjj  is  bound  and  stress 
singularity  no  longer  exists  at  the  tip. 


IV.  One-dimensional  Approximate  Analysis  of  Non-local  Tensile 
Stress  in  Front  of  Crack 


In  order  to  obtain  a  clear  qualitative  physical  picture  of 
the  non-local  stress  field  at  the  tip,  in  this  section  we  will 


atteopt  to  simplify  the  equation  to  one-dimension,  i.e.,  to  use 


the  following 


/(*) 


(23) 


to  calculate  the  tensile  stress  in  front  of  the  crack  (for 
convenience,  subscripts  are  omitted).  This  is  equivalent  to 
assuming  that  the  stress  along  x2-direction  remains  unchanged  at 
least  in  the  atomic  scale  and  non-local  action  is  only  active 
along  the  crack. 

From  equations  (20)  and  (23)  it  is  easy  to  find 

o  when  0<x<a-b 


(!+«?/«-« 
2+*  \  b 


when  a-b<x<a 


^  "  I  2  +  ■  [(~T~  +  0^  ”  2(~~»  ")^]  *  when  a<x<a+b 


a + »y  r 

(  *—  •  _t_  ,Yfsr  l(*  +  1 

r*-*.  j\f6i 

2  +  »  l 

{  b  l)  \ 

when  x>a+b 


where 
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Figures  4  and  5  show  the  tensile  stress  distribution  with 
different  internal  structure  scales. 


Figure  4.  Tensile  Stress  Distribution  in  Front  of  Crack  when 
a/b=50.  Dotted  line  represents  the  classical 

solution  when  n*0.2,  t  /o  =0.5. 

00  y 


Figure  5.  Tensile  Stress  Distribution  in  Front  of  Crack  when 
a/b  *200 . 


It  is  obvious  that  as  long  as  b  is  not  zero,  the  stress  in 


front  of  the  crack  is  always  finite.  The  maximum  stress  occurs 
at  r»b/21+n-1.  It  is 

in  -  A.  [-0  (25) 

where 

.  (1  +  «)W0)  r*  091  — 

+  ■  +  0J1 

is  a  constant  related  to  the  hardening  index,  as  shown  in  Figure 
6.  Based  on  this,  the  stress  concentration  factor  is 

(26) 

It  is  not-  only  related  to  the  relative  scale  of  crack  to  internal 
structure,  but  also  to  the  hardening  index  as  well  as  to  the 
relative  magnitude  of  external  load  to  yield  stress.  This  is  an 
apparent  result  due  to  non-linearity  of  plasticity. 


Figure  6.  Coefficient  A  vs.  Hardening  Index  n. 


Since  the  stress  near  the  tip  is  finite,  naturally  we  can 

establish  a  physical  criterion  for  cracking:  there  is  a  critical 

stress  t„ ,  which  is  a  constant  for  a  material,  when  t  =t  , 
c  ’  '  max  c  ’ 

crack  will  begin  to  propagate. 

Equation  (25)  can  also  be  written  as 


-  As, 


where  E  is  the  Young's  modulus.  Because  n,  A  ,  b,  o  ,  and  E  are 

n  y 

material  constants,  therefore,  the  criterion  in  classical  / 

elastic-plastic  fracture  mechanics  can  be  directly  obtained  based 
on  this  fracture  criterion.  Thus,  the  contradiction  in  classical 
fracture  mechanics  as  pointed  out  in  the  introduction  can  be 
bypassed.  The  macroscopic  model  of  fracture  mechanics  is 
unified  with  the  microscopic  physical  mechanism. 

It  is  worthwhile  to  note  the  physical  significance  of 
critical  stress  tc.  The  quasi-embrittlement  fracture  and  ductile 
fracture  of  alloys-  are  far  more  complicated  than  the  fracture  of 
a  single  crystal.  On  one  hand,  because  of  the  presence  of 
intergranular  boundary,  alloying  elements  and  various 
metallographic  structures,  the  microscopic  structure  is  non- 
uniform.  On  the  other  hand,  as  plastic  deformation  progresses, 
the  microscopic  structure  also  continues  to  change.  Despite  so, 
we  can  still  anticipate  that  the  critical  stress  t  corresponds 
to  the  case  that  the  cohesive  strength  between  atoms  at  the 
leading  edge  of  the  crack  has  the  same  order  of  magnitude  as  the 
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cohesive  strength  of  the  lattice  of  the  base  metal  (which  is 
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around  0.1E  for  most  metals)  after  undergoing  a  severe  plastic 
deformation.  Table  1  shows  the  tc/E  values  for  several  steel 
materials  calculated  from  available  experimental  data.  The 
results  are  in  good  agreement  with  the  above  analysis.  From  the 
data  one  can  see  that  although  materials  of  similar  composition 
may  have  large  variations  in  yield  strength  a ,  yet  their  tc 
values  are  relatively  stable  when  the  fracture  mode  is  identical 
The  numerical  value  is  slightly  less  than  the  cohesive  strength 
of  a  perfect  lattice.  This  reduction  is  not  surprising  if  we 
consider  factors  such  as  local  stress  concentration  due  to  a 
large  number  of  lattice  defects  and  dislocations  near  the  tip 
prior  to  failure.  Obviously,  the  microscopic  structure  of  the 
material  is  different  due  to  various  heat  treatment  conditions. 
Therefore,  the  effect  of  the  above  factors  is  also  different. 
Thus,  the  value  of  t  varies  slightly.  In  Table  1,  the  data  on 
two  Fe-20%  Co- 15%  Cr-5%Mo  alloys  reflect  this  variation.  The 
data  are  arranged  i-n  descending  order  based  on  effective  time.  A 
it  increases,  the  residue  Austinite  content  decreases  and  the 
precipitate  phase  increases.  Thus,  the  atomic  binding  force 
decreases  and  the  local  stress  concentration  effect  is 
strengthened,  t  will  decrease  slightly. 


’•  tor?”ng'sCSid“«Bl”  Several 

Materials 


m  # 

(MN/nO 

9 

L 

CMN/«) 

W* 

636 

0.17 

0.086 

0.093 

GCr  15 

431 

0.21 

0.046 

0.104 

[13] 

* 

411 

0.21  . 

0.030 

0.091 

34  CrNiMo 

1225 

0.10 

0.032 

0.060 

... 

.  6  j 

1450 

0.085 

0.038 

0.059 

980 

0.12 

0.119 

0.073 

[6] 

30SiMnCrMo 

1156 

0.11 

0.117 

0.073 

[1-1 

1156 

0.10 

0.119 

0.064 

• 

1390 

0.066 

Q.C413 

0.046 

F«-20%Co-I596Cr- 

1380 

0.058 

0.C247 

0.045 

5%Mo 

1750 

0.051 

0.0203 

0.045 

[15] 

3.  C*«tt 

2110 

0.0J4 

0.0057 

0.042 

2300 

0.021 

o.eas 

0.040 

1590 

0.076 

0.0456 

0.058 

Fe-20%Co-1596Cr- 

5%Mo 

1510 

0.070 

0.0215 

0.U58 

[15] 

a-.  CXA.S40 

2140 

0.050 

0.0051 

0.050 

2360 

0.033 

0.0018 

0.045 

1.  type  of  steel 

2.  source  of  data 

3.  high  purity 

4.  high  purity,  de-oxygenated 

V.  Discussion  and  Conclusions 
1.  In  spite  of  the  fact  that  the  structure  is  only  given 
based  on  plastic  deformation  in  this  paper,  this  basic  concept  is 
also  applicable  to  other  types  of  plastic  structure  relations. 
Regardless  which  structural  assumption  is  adopted,  we  must  use  I 
strain  space  to  describe  the  problem.  In  this  space,  the 


classical  stress  vs.  elastic  strain  relation  can  be  used  to 
transform  stress  into  elastic  strain.  This  type  of  treatment 
also  ensures  that  the  present  model  is  consistent  with  the 
classical  model  when  the  internal  scale  in  the  non-local  modulus 
approaches  zero. 

2.  Different  from  the  non-local  plastic  theory  recently 
r  i6] 

introduced  by  Eringen  ,  this  model  (does  not  take  the  non¬ 
local  characteristic  of  the  plastic  structure  into  account.  This 
is  because  the  physical  characteristic  scale  relating  to  plastic 
deformation  is  far  larger  than  the  characteristic  scale 
considered  in  non-local  elasticity  (atomic  scale).  Relative  to 
atomic  scale,  plastic  deformation  is  highly  non-uniform.  From 
another  angle  one  can  say  that  the  plastic  deformation  of  a  point 
is  somewhat  random  in  nature.  Hence,  when  the  accuracy  is  on  the 
atomic  scale  we  should  understand  that  the  continuity  of  plastic 
strain  is  described  as  the  statistical  average  of  all  probable 
deformations  in  the  neighborhood  of  the  point  under 
investigation.  The  physical  connection  between  plastic  strain 
and  stress,  however,  is  expressed  as  the  internal  stress  caused 
by  the  growth  of  dislocations.  Statistically,  it  is  equal  to 
zero.  Therefore,  macroscopically  plastic  strain  is  not  directly 
related  to  stress.  However,  microscopically,  this  internal 
stress  causes  the  average  cohesive  strength  to  drop. 
Mathematically,  the  model  introduced  in  this  work  decouples  the 
non-local  (elastic)  and  non-linear  (plastic)  parts  of  the  basic 
equations.  Consequently,  it  is  more  convenient  to  use  in 


practical  problems. 

3.  Strictly  speaking,  in  the  fracture  region  near  the  tip, 
geometrical  non-linearity  caused  by  a  large  strain  cannot  be 
neglected.  J  integral  is  only  applicable  outside  this  region. 
Thus,  if  the  entire  region  is  treated  based  on  the  small  strain 
linear  theory,  it  will  read  to  some  error.  Despite  so,  the  non¬ 
local  elastic  plastic  theory  recommended  in  this  work  can  unify 
the  macroscopic  model  of  fracture  mechanics  with  the  microscopic 
physical  mechanism.  Furthermore,  the  critical  fracture  stress 
thus  obtained  is  on  the  same  order  of  magnitude  as  the  cohesive 
stress  of  the  lattice.  This  will  benefit  the  understanding  of 
the  physical  mechanism  of  the  material  fracture  process. 

4.  In  metal  physics,  the  effect  of  microstructure  factors 
such  as  two-phase  particle  volume  integral,  particle  size  and 
impurity  spacing  on  the  fracture  resistance  is  discussed  in 
detail.  In  the  present  model,  these  factors  affect  the  fracture 
resistance  of  materials  through  Oy  and  n,  i.e.,*the  change  of 
plastic  deformation  in  front  of  the  crack. 

5.  The  fracture  of  a  metal  involves  three  different  scales: 
macroscopic,  fine  detail  and  microscopic.  In  this  work,  an 
attempt  was  made  to  connect  microscopic  and  macroscopic 
viewpoints.  From  the  fine  details,  many  physical  factors  are  yet 
to  be  considered.  It  is  necessary  to  understand  the  dependence 
of  internal  cohesive  stress  upon  fine  factors  such  as  alloy 
composition,  metallograph  and  microscopic  structure  as  well  as 
upon  macroscopic  parameters  such  as  hardening  index  and 


destructive  strain  which  describe  the  extent  o£  plastic 
deformation.  A  great  deal  of  detailed  theoretical  and  experimental 
work  is  required. 
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Abstract 

A  model  of  nonloeal  elastic-plastic  oontinuum  is  proposed.  The  stress  and  the  elastic 
strain  are  related  by  a  nonlocal  linear  relation,  and  the  plastic  strain  is  dependent  on 
the  history  of  total  strain.  For  plastic  deformation  theory,  it  is  assumed  that  the  plas¬ 
tic  strain  tensor  is  proportional  to  the  total  strain  deviation  tensor  and  the  ratio  is  a 
sealer  function  of  the  effective  total  strain. 

This  model  is  used  to  analyse  the  stress  field  at  the  tip  of  a  crack  in  a  power-law 
hardening  material  under  plane  strain  condition.  $ased  on  the  results  of  HRR  asymp¬ 
totic  solution  in  classical  fracture  mechanics,  the  distribution  of  tensile  stress  on  the  li¬ 
ne  directly  ahead  of  a  crack  tip  and  the  expression  for  maximum  tensile  stress  are  ca¬ 
lculated  under  one-dimensional  simplification.  It  is  shown  that  the  J.  criterion  may  be 
obtained  from  the  maximum  tensile  stress  criterion  of  the  nonlocal  theory.  Available  ex¬ 
perimental  data  for  steels  are  used  to  calculate  the  maximum  tensile  stress  of  free  tut  e 
initiation  at  the  crack  tip.  which  is  found  to  be  done  to  the  theoretical  cohesive  strength. 
The  results  obtained  an  useful  for  understanding  the  fracture  process  and  mechanism 
of  materials. 


Study  of  Plane  Stress  with  Elastic  Plastic  Mixed  Mode  Fracture  ' 
Xu  Jilin,  Xue  Yinian  and  Han  Jinhu 
(Institute  of  Mechanics,  Academia  Sinica) 

Abstract 

In  this  work,  the  deformation  field  around  the  crack  of  a 
thin  aluminum  plate  with  center  cracks  of  various  angles  of 
inclination  under  tensile  load  was  measured  using  the  direct  laser 
speckle  method  and  Moire’s  method.  In  addition,  the  relation 
between  loading  and  crack  propagation  in  the  steady-state  crack 
growth  process  was  measured.  Moreover,  finite  element  analysis 
based  on  large  elastic-plastic  deformation  equation  was  carried 
out  to  obtain  the  stress-strain  distribution  around  the  crack. 
Calculated  results  are  found  in  good  agreement  with  experimental 
data.  Furthermore,  the  results  were  discussed. 

I.  Introduction 

In  a  thin  wall  structure  of  a  malleable  material  with  a 
crack,  before  the  crack  began  to  propagate,  the  crack  tip  has 
already  formed  a  large  plastic  region.  After  the  crack  is 
Initiated,  there  is  a  steady-state  growth  process.  If  the  linear 
elastic  fracture  theory  is  still  applied,  the  load  capability  of 
the  structure  will  be  under-estimated.  Therefore,  there  is  a 
need  to  establish  an  elastic-plastic  fracture  theory  to  conduct 
this  investigation.  In  recent  years,  progress  has  been  made  in 
the  study  of  elastic-plastic  mode  I  fracture  theory.  New  methods 


and  parameters  are  constantly  being  Introduced.  In  references 

[1,2],  many  authors  used  different  fracture  criteria  to  study  the 

Initiation,  steady-state  propagation  and  ductility  instability 

over  a  wide  range  of  yield  conditions.  In  addition,  experimental 

data  and  finite  element  numerical  analysis  were  combined  to 

calculate  the  stress-strain  field  near  the  crack  tip. 

Furthermore,  various  fracture  parameters  were  analyzed.  Among 

them,  Shih  et  al  conducted  a  great  deal  of  experimental  work  and 

numerical  analysis  under  plane  strain  for  large  range  yield. 

After  sorting  various  fracture  parameters,  it  was  recommended 

that  Jlc  and  6X  be  used  to  express  the  initiation  of  a  crack, 

and  £he  tear  modulus  Tj  and  Tfi  be  used  to  express  the  growth  of 

the  crack.  Kanninen  believed  that  J  and  CTOA  are  the  most 

effective  parameters  to  estimate  the  initiation,  steady-state 

growth  and  ductility  instability  of  a  malleable  material.  In  the 

[3] 

study  of  the  plane  stress  fracture  problem,  Feddersen  used 
aluminum  alloy  plates  with  center  cracks  to  conduct  a  large 
number  of  experiments.  A  useful  analytical  method  for 
engineering  design  was  provided.  However,  the  steady-state  crack 
propagation  process  was  only  quantitatively  described.  We 
measured  the  relation  between  steady-state  crack  growth  and  load 
increment  by  using  thin  aluminum  alloy  plates  with  center  cracks 

[A] 

under  tensile  load  .  Moreover,  based  on  the  elastic  plastic 
fracture  model  of  plane  stress  —  the  shrink  neck  band  model^^, 
the  elastic-plastic  deformation  finite  element  method  was  used  to 
calculate  the  steady-state  growth  of  Mode  I  crack^^.  The 
calculated  results  and  experimental  data^^  are  in  good 


agreement.  There  is  little  published  work  on  the  study  on  mixed 
mode  fracture  for  large  range  yield.  Ueda  et  al^^  studied  the 
initiation  of  Mode  I  and  Mode  II  cracks  under  large  range  yield 
conditions.  They  conducted  some  low  temperature  embrittlement 
tests  with  soft  steel  cross  specimens  with  inclined  cracks  under 
the  tensile  stress  in  two  axes.  In  the  meantime,  they  conducted 
elastic-plastic  finite  element  calculation  for  plane  stress. 

They  used  the  COD  concept  for  Mode  I  fracture  to  analyze  the 
initiation  of  a  mixed  mode  crack. 

In  this  work,  thin  aluminum  alloy  plates  with  center  cracks 
at  various  angles  of  inclination  under  tensile  load  were  used  in 
the  experimental  study  in  plane  stress 
conditions.  Both  Ithe  laser  speckle  method  and  Moire’s 
used.  In  addition,  the  finite  element  method  with  elastic-plastic 
deformation  was  used  to  perform  numerical  analysis.  The  results 
were  compared  to  experimental  data  and  discussed. 

II.  Experimental 

This  experiment  involves  20  specimens.  They  are  thin  plates 
with  center  cracks  made  of  LY12-CZ  and  LY12-CS  aluminum  alloys. 

Manuscript  received  on  October  15,  1983.  The  paper  was 
presented  at  the  1983  Beijing  International  Fracture  Mechanics 
Academic  Meeting. 

The  characteristic  data  of  the  materials  are  shown  in  Table  1.  / 

These  numbers  were  obtained  in  the  tensile  test  of  the  raw 
materials.  In  the  table,  a  and  n  are  parameters  in  the 
expression  e-  a/E  in  which  an  index  curve  is  used 


to  approximate  the  curve  of  the  material.  The  geometric 
dimensions  of  the  specimen  are  shown  in  Figure  1 .  The  angle  of 
inclination  between  the  crack  line  and  the  load  B=90°,  60°,  45°, 
and  30°.  The  initial  crack  length  is  2aD.  It  ranges  from  4.65mm 
to  42.7mm.  Their  projections  in  the  direction  perpendicular  to 
the  load  range  from  4.65mm  to  23.8mm. 

Table  1 


1 .  material 

2.  elastic  modulus 

3.  yield  strength 

4.  limiting  strength 

5.  elongation 


Figure  1 

The  center  crack  In  the  sample  is  cut  by  a  0.12mm  diameter 
molybdenum  and  made  into  a  fatigue  crack.  Fatigue  cracks  were 
not  prepared  on  six  specimens  in  order  to  study  the  effect  of 
prepared  fatigue  crack  on  crack  propagation. 

Because  the  specimen  is  very  thin,  there  is  not  much 
difference  between  surface  and  internal  crack  growth.  In  the 
experiment,  an  80X  microscope  was  used  to  read  the  crack  growth 
increment  Aa.  The  minimum  scale  on  the  microscope  is  0.01mm. 


The  Initiation  and  steady-state  growth  can  be  clearly  observed 
through  the  microscope.  The  relation  between  load  o-and  crack 
growth  (a-a, )  is  thus  determined. 

In  order  to  measure  the  strain  distribution  around  the  crack 

during  the  entire  loading  process,  a  laser  speckle  method  is  used 

to  directly  measure  the  deformation  on  one  surface  when  the  load 

is  relatively  low.  On  the  other  surface,  Moire's  method  is  used 

[81 

to  determine  the  deformation  when  the  load  is  relatively  high  . 
The  surface  on  which  the  direct  laser  speckle  method  is  used  is 
polished  by  a  wheel  to  Improve  the  reflectivity  of  the  surface  in 
order  to  obtain  high  quality  double  exposure  film  for  the 
analysis  of  the  entire  field.  The  analytical  fringe  pattern  of 
the  entire  field  is  photographed  through  a  2mm  diameter  filter 
hole.  The  hole  is  opened  up  high  to  reach  a  sensitivity 
corresonding  to  a  grating  line  density  of  467  lines  per 
millimeter.  Hence,  it  is  possible  to  measure  the  deformation 
field  when  the  load  is  relatively  low.  When  the  load  is 
relatively  high,  a  simple  Moire's  method  can  be  used  to  measure 
on  the  other  surface  even  when  deformation  is  large.  The  grating 
line  density  chosen  is  40  lines  per  milliliter  and  the  specimen 
grating  is  an  orthogonal  grating.  The  analyzing  grating  is  a 
unidirectional  grating  of  the  same  density.  Based  on  the 
experimental  data  shown  in  Figure  8,  the  results  of  the  speckle 
method  and  those  of  Moire's  method  are  in  good  agreement. 

In  order  to  investigate  whether  the  elastic-plastic  Mode  I 
fracture  model  of  plane  stress^^  and  crack  growth  criterion"^ 


•!'.  'I'.r.'vr^ 1  -  w.r.^'"rf 


can  be  extended  to  mixed  mode  fractures ,  we  are  concerned  about 

the  strain  distribution  on  the  ductile  belt  of  the  Inclined 

center  crack  specimens  (extension  from  the  two  apices  of  the 

crack  along  the  direction  perpendicular  to  the  tensile  direction  is 

called  the  ductile  belt).  Moire's  method  was  used  to  measure  the 

strain  e  in  the  x-direction  and  strain  e  in  the  y-direction. 
x  y 

Results  show  that  e  is  much  smaller  than  e  . 

x  y 

HI.  Finite  Element  Analysis  / 497 

In  this  work,  the  Euler  finite  element  formula  with  plane 

stress  deformation  was  used  to  calculate  and  analyze  aluminum 

alloy  plates  with  center  cracks  under  load  at  angles  of 

inclination  6=60°,  45°  and  30°  to  determine  the  stress  o  ,  o  ,t 

x  y  xy 

in  the  loading  process,  the  stress  o  .oa  and  t  -  around  the  crack 

r  o  n? 

top,  the  strain  e  .  c  ,  e_,  y  and  the  distribution  of  the 

x  y  z  xy 

plastic  region.  In  addition,  the  open  displacement  of  the  crack 
surface  and  the  strain  on  the  ductile  belt  were  calculated  and 
compared  to  the  experimental  result  (see  Figures  8  and  9). 

With  regard  to  the  steady-state  crack  growth  on  thin 
aluminum  alloy  plate  specimens  with  center  cracks  at  0=90°,  we 
had  calculated  the  relation  between  crack  growth  and  load 
increment,  as  well  as  the  instability  load,  using  the  finite  element 
method  with  plastic  deformation  equation  based  on  the  shrink  neck 
band  model  of  plane  stress  elastic-plastic  fracture^^.  The 
criterion  of  crack  propagation  was  that  the  relative  elongation 
of  the  shrink  neck  region  at  the  crack  tip  reached  the  same  level 
as  the  elongation  coefficient  of  the  material The  results  are 
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in  good  agreement  with  the  experimental  data. 

The  speed  equilibrium  equation  in  the  form  of  pseudo-work 


equation  is 


- j  <r„S(2DnDi,  —  dV 

—  |  f.SVjds  +  iiStjdV 


(1) 


where  is  the  Jaumann  speed  of  Kirchhoff  stress  tensor  aij 

and  represent  the  true  stress  tensor  and  deformation  speed 

tensor,  respectively,  f^  and  are  the  instantaneous  unit 

* 

surface  force  and  unit  volume  force,  and  are  related  by 

the  constitutive  matrix  of  elastic  deformation  increment  [C]  • 

(2) 

From  equation  (1),  we  get  the  Euler  finite  element  speed 
equilibrium  equation 

[kh<m  -  jy[<v]WF  +  [mm}d-  (3) 


mmm 

Baa 


Figure  2 
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4}  is  the  nodal  point  speed  array.  The  rigidity  matrix  is 

[K]  -  (  {[BHC1IB]  +  -  2[Bkiyffii[BKt])dV  (4) 

In  the  plane  stress  case,  let  us  assume  that  a  =a  =o  =0.  The 
r  1  x  xx  yz 

instantaneous  element  thickness  t=t0  eEx  . 

When  an  increment  tangent  rigidity  method  is  used  to  solve 
the  problem,  the  load  increment  coefficient  is  determined  by 
using  the  "Shantian"  method  to  allow  each  element  to  enter  a 
yield  state  in  order  to  reduce  the  error  introduced  due  to  the 
discontinuity  of  the  tangential  modulus  of  the  curve  of  the 
material  at  the  elastic-plastic  turning  point.  Moreover,  it 
corrects  for  an  unbalanced  load. 

-  The  unit  mesh  for  an  inclined,  center-cracked  specimen  is 
shown  in  Figure  2.  The  center-cracked  unit  lattice  whose  ratio 
of  unit  dimension  at  the  crack  tip  to  the  crack  length  is 

1/40-1/100,  0=90 

is  shown  in  reference  [6], 

IV.  Results  and  Discussion 

1.  Fracture  Surface  Topography 

From  the  fracture  surfaces  of  all  failed  specimens  one  can 
see  that  there  is  a  very  small  triangular  flat  cross-section  in 
the  neighborhood  of  the  tip.  It  is  rapidly  turned  into  a  45° 
shear  fracture  with  respect  to  the  surface.  It  belongs  to  a 
plane  stress  fracture. 

2.  Fracture  Stress 

Table  2  shows  the  non-dimensional  experimental  data 
including  the  uniform  tensile  stress  o.  and  average  stress 


°l(net)  on  t*ie  ductile  belt  at  initiation  and  the  uniform  tensile 
stress  oc  and  average  stress  °c(net)  on  the  ductile  belt  during 
unstable  growth.  These  results  show  that  the  average  stress  on 
the  ductile  belt  prior  to  the  unstable  growth  of  the  crack  has 
exceeded  the  yield  of  the  material.  Therefore,  they  most  failed 
under  the  condition  of  large  range  yield. 


Figure  3a 

(1)  specimen  28 

(2)  specimen  24 

(3)  specimen  19 
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Figure  3b 

(1)  specimen  7 


During  initiation  and  unstable  growth,  the  uniform  tensile 


stress  and  oc  increase  with  increasing  W/b  (ratio  of  ductile 
belt  width  to  plate  width).  For  different  angle  of  inclination 
ft  as  long  as  the  ductile  belt  width  to  plate  width  ratio  W/b  is 
close,  i.e. ,  when  the  projection  of  the  crack  along  the  direction 
perpendicular  to  the  load  2a<,(a2asinft)  is  close,  the  fracture 
stress  a  is  almost  identical.  Furthermore,  in  plane  stress 
conditions,  the  average  tensile  stress  during  unstable  growth  a 

v 

is  basically  independent  of  the  plate  thickness. 

3.  Pre-f abricated  Fatigue  Crack 

From  the  experimental  data  obtained  with  six  specimens 
without  pre-fabricated  fatigue  cracks  shown  in  Table  2  one  can 
see  that  their  crack  initiation  load  is  very  close  to  their 
instability  load  aa.  There  is  a  lack  of  an  obvious  steady-state 
crack  growth  stage.  The  instability  load  of  specimens  with  pre¬ 
fabricated  fatigue  cracks,  however,  is  apparently  higher  than  the 
initiation  load,  furthermore,  there  is  an  apparent  steady-state 
crack  growth  process  as  shown  in  Figure  3. 

4.  Crack  Propagation  Process 

Figure  3a  shows  the  results  of  the  steady-state  crack 
propagation  process  of  three  specimens  (No.  19,  24  and  28)  at  3 

=60° ,  45°  and  30° ,  including  the  relation  between  uniform  tensile 

/  50C 

stress  o  with  crack  growth  increment  Aa.  Figure  3b  shows  the  o- 
Aa  relations  measured  experimentally  as  well  as  calculated  using 
the  finite  element  method  based  on  Mode  I  fracture  (specimen  7) 
at  8=90°^^.  The  calculated  a  is  in  good  agreement  with  the 


experimental  data.  The  error  is  less  than  5%. 
5.  Direction  of  Crack  Propagation 


For  center-cracked  thin  plate  specimens  with  different 
angles  of  inclination  B  under  tensile  load,  the  crack  propagates 
essentially  in  the  same  manner  as  that  of  a  Mode  I  crack.  It 
grows  steadily  along  a  direction  perpendicular  to  the  tensile 
direction  until  reaching  an  unstable  state.  Figure  4  shows  the 
photographs  of  three  failed  specimens  at  8=60°,  45°  and  30°.  Th 
angle  between  the  direction  of  crack  propagation  and  the 
direction  perpendicular  to  the  load  line  (x-axis)  is  less  than 
10°. 


■  4 


Figure  4 
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Figure  5.  Specimens  19,  24  and  28 

(Values  of  o/o  are  0.728,  0.711  and  0.712, 
respectively)  y 

Figure  5  shows  the  tensile  stress  distribution,  aQ,  around 
the  crack  tip  near  the  initiation  load  as  calculated  by  the 
finite  element  method.  In  the  element  where  the  angle  BE  with 
the  y  axis  is  close  to  90°,  oQ  reaches  its  maximum  value.  This 

w 

•  X 

coincides  with  experimental  observation  of  crack  propagation 

perpendicular  to  the  load  direction. 

In  addition,  based  on  the  displacement  field  calculated,  the 

relative  displacement  of  the  point  approximately  0.4mm  away  from 

the  crack  top  is  named  the  crack  top  opening  vector  c5Dr  (the 

numerical  value  of  c5d  is  /  u* +  u*  where  u  and  u  are  relative 

r  x  y  x  y  _ 

displacements  in  x  and  y  direction).  The  direction  of  C0Dr 
obtained  from  finite  element  calculation  is  near  the  load  line 
direction.  For  specimens  with  B=60° ,  45°  and  30°,  the  angle 
between  C0Dr  and  the  load  line  is  approximately  4°-10° .  This 


fact  indicates  that  the  estimated  direction  of  crack  propagation 
is  essentially  consistent  with  the  experimental  results. 


6.  Plastic  Region 


Figure  6.  Specimens  19,  24  and  28 

(o/o  values  are  0.728,  0.711  and  0.712, 
respectively) 


Figure  7.  Specimens  19,  24  and  28 

Figure  6  shows  the  plastic  region  distribution  as  calculated 
by  the  finite  element  method.  Before  the  crack' begins  to 
propagate,  the  plastic  region  is  already  not  quite  small.  The 
plastic  regions  in  specimens  with  various  angles  of  inclination  3 
are  similar  in  shape.  The  plastic  region  width  on  the  ductile 
belt,  Xp,  is  given  in  Figure  7.  The  calculated  values  obtained 
from  finite  element  method  agree  well  with  the  experimental  data. 

7.  Strain  Field 


The  y-direction  strain,  e^,  on  the  ductile  belt  is  shown  in 
Figures  8  and  9.  Results  obtained  by  the  laser  speckle  method  and 
Moire's  method  agreed  very  well.  In  the  figures,  the 
experimental  values  are  generally  lower  than  the  calculated 
results,  primarily  because  the  strain  measured  is  an  average 
between  fringes  and  usually  the  spacing  is  larger  than  the 
dimension  of  the  mesh  used  in  a  finite  element  method. 

Both  experimental  and  theoretical  results  show  that  the  x- 
direction  strain,  ex,  is  much  smaller  than  the  strain  in  the 
tensile  line  y-direction,  e^,  for  specimens  with  various  angles 
of  inclination.  The  computation  also  shows  that  shear  strain  on 
the  ductile  belt,  y  ,  is  also  much  smaller  than  e  .  Therefore,  ey 
is  the  primary  strain  on  the  ductile  belt. 

In  summary,  under  plane  stress  conditions,  for  center 
cracked  specimens,  the  fracture  stress  is  primarily  determined  by 
the  ductile  belt  width  to  plate  width  ratio,  W/b.  This  means 
that  the  fracture  stress  is  determined  by  the  projection  of  the 
crack  in  a  direction  perpendicular  to  the  load  line,  2a0  . 
Furthermore,  it  is  basically  independent  of  the  plate  thickness. 

In  center  cracked  specimens  of  various  angles  of  inclination  0 
under  uniform  tensile  load,  the  crack  basically  propagates 
perpendicular  to  the  load  line.  The  angle  e  between  the 
direction  of  cracking  and  the  perpendicular  direction,  is  less 
than  10°.  After  the  crack  is  initiated,  within  a  small  range 
where  the  crack  growth  increment  is  less  than  the  plate 
thickness,  0  rapidly  approaches  0°.  In  addition,  e^,  the  strain 
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in  the  y-direction,  is  the  primary  strain  on  the  ductile  belt. 


These  facts  indicate  that  it  is  possible  to  adopt  the  Mode  I 
elastic-plastic  fracture  raodel^^  and  the  criterion  of  crack 
propagation  in  which  the  relative  elongation  reaches  the 
elongation  of  the  material^^  to  solve  the  plane  stress  mixed 
mode  fracture  problems . 
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AN  E1AST1C-PLAST1C  MIXED  MODE  FRACTURE 
-  INVESTIGATION  FOR  PLANE  STRESS 

Xu  Jilin,  Xne  Yinian,  Han  Jinhu 
(Institute  of  Ueduoua,  Acodenrit  Shoes) 

Abstract 

In  thia  paper  the  stable  crack  growth  processes  in  aluminium  alloy  sheet  apeimens 
with  flat  or  inclined  central  crack,  subjected  to  uniform  tensile  load  have  been  inves¬ 
tigated.  The  relation  between  the  tensile  load  and  the  amount  of  crack  extension  was 
obtained.  The  deformation  field  was  measured  by  using  laser  speckle  mehtod  and 
Moire’  method.  Finite  element  analysis  baaed  on  large  elastic-plastic  deformation  equa¬ 
tion  has  also  been  carried  out.  The  criterion  proposed  to  predict  crack  growth  under 
mode  I  condition  is  that1*1  the  tensile  strain  at  the  crack  tip  reaches  the  maximum 
elongation  of  material.  The  calculated  resulte  are  in  good  agreement  with  the  experi¬ 
ment  data. 
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An  Analytical  Solution  of  Dynamic  Response  for  Ideal  Rigid 

Plastic  Timoshenko  Beam 
Jin  Quanlin 

(Research  Institute  of  Mechanical  and  Electrical  Technology, 
Ministry  of  Machine  Building  Industry) 

Abstract 

In  this  paper,  an  analytical  solution  for  a  fixed-ended 
Timoshenko  beam  under  uniform  dynamic  load  is  given  by  using 
different  discontinuity  conditions  for  various  states  of  motion 
at  the  rigid  plastic  interface.  This  solution  is  applicable  to 
any  non-reversing  load  which  varies  with  time.  At  the  end,  the 
effect  of  the  rotational  inertia  on  the  dynamic  response  of  the 
beam  is  discussed. 


I.  Introduction 

In  recent  years  quite  a  few  people  have  studied  the 
influence  of  shear  effect  and  rotational  inertia  on  the  plastic 
dynamic  response  of  the  beam1  ”  .  Jones1  ”  summarized  these 

results.  However,  their  studies  are  limited  to  dealing  with  the 
initial  velocity  problem  which  is  equivalent  to  the  immediate 
unloading  after  instantaneous  loading.  It  is  very  important  to 
use  the  discontinuity  condition  accurately.  To  this  end, 
reference  Cl]  was  the  first  paper  to  discuss  it  theoretically. 

It  was  pointed  out  that  different  discontinuity  conditions  should 
be  used  for  different  directions  of  motion  on  the  discontinuity 
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plane.  In  this  paper,  a  fixed-ended  Timoshenko  beam  under  a 
uniform  load  is  used  as  an  example  to  solve  the  problem  states  of 
motion  at  a  rigid  plastic  interface.  Furthermore,  the 
transformation  condition  from  one  state  to  another  is  given. 
Therefore,  the  solution  given  in  this  work  is  applicable  to  any 
non-reversing  time-varying  load. 


■  1  ■  2 


Figure  1  Figure  2 

II.  Basic  Equations 

Let  us  consider  a  fixed-ended  (neglecting  axial  force), 
equi-sectional  ideal  plastic  beam  which  is  1  in  length  and  is 
under  a  uniform  dynamic  load  q(t)  (see  Figure  1).  Its  mass  per 
unit  length  is  m,  and  the  rotational  inertia  is  r.  The  bending 
moment  and  shear  force  on  the  cross-section  of  the  beam  are  M  and 
Q,  respectively.  (The  positive  directions  are  shown  in  Figure 
2).  The  defelction  of  the  beam  is  y.  The  angle  of  inclination 
caused  by  the  bending  deformation  of  the  beam  is  4>.  x  is  used  to 
represent  the  abscissa  of  the  beam  cross-section.  The  equation 


of  motion  of  the  beam  is 


Q'  —  m9  ~  f»M'  ~  Q  ~  r*r24>.  (1) 

where  the  dots  on  top  of  the  letters  and  the  quotation  mark  on 
the  upper  right  corners  of  the  letters  are  the  partial 
derivatives  of  time  t  and  x.  The  total  angle  inclination  of  the 
beam  is  y,=4>+y,  where  y  is  the  rotation  angle  caused  by  shear 
deformation.  Furthermore,  it  is  specified  that  the  rotation  of  <|> 
and  y  in  the  positive  direction 
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makes  y’>0. 


Let  us  introduce  the  following  non-dimensional  quantities: 
£  -  M/M„  Q  -  QIQ ..  *  -  fj1/ 2M.»  P  “  £>»//«„  (y)  , 


ml2  T  _  ml 3  ,  t 

nT.”  * 


fnl*  t  -  ml*  T  —  x  f  —  ^  M  * 

ru,2’  '  ru.T’  ’  T’c  i’ ‘  T ■ 


where  M*  and  Q0  are  the  limiting  bending  moment  of  pure  bending 
yield  of  the  beam  cross-section  and  the  limiting  shear  force  of 
pure  shear  yield,  respectively.  A  is  the  x-coordinate  at  the 
rigid  plastic  interface.  K=<j»’  is  the  curvature  of  bending  of  the 
beam.  T  is  a  unit  time.  In  the  following,  for  convenience,  the 
bars  on  top  of  M,  Q,  t,  -x,  4.,  k  and  y  are  omitted.  They  are 
expressed  as  non-dimensional  quantities.  Equation  (1)  is  re¬ 
written  as: 

M'-pQ-ai  (2) 


Let  us  choose  a  square  yield  plane  as  shown  in  Figure  3.  In  the 
plastic  region,  equation  (2)  is  solved  by  using  the  yield 
condition  and  the  relevant  orthogonal  flow  method.  When  stress 
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points  are  located  on  the  side  AB  or  CD  of  the  yield  plane, 

•  —  CXO  A  -7==  +  C,(/)  ch  — —=  +  2^.  (3  ) 

.Ve  V®  ' 

When  the  stress  points  are  located  on  the  side  AC  or  BD  on  the 
yield  plane, 

y»2u  (4) 

In  the  rigid  region,  we  have 

y»C,  (t)x+C4  (t)  (5) 

In  order  to  determine  (^  ,  Ca  ,  Ca  and  C4  in  these  equations,  it  is 
necessary  to  know  the  boundary  condition  of  the  beam  and  the 
following  discontinuity  condition^ ^  at  the  rigid  plastic 
interface: 

.When. the  motion  of  the  rigid  plastic  interface  leads  to  the 
expansion  of  the  plastic  region, 

[y]-[VJ=0  and  [Q]*[M]*0  (6) 

When  the  motion  of  the  rigid  plastic  interface  leads  to  the 
contraction  of  the  plastic  region, 

[y]«U3=0  and  CQ3=[M]=0  .  (7) 


[z]=z+-z“  which  is  the  discontinuity  value  of  a  physical  quantity 
z  at  the  rigid  plastic  interface.  z+  and  z"  are  the  values  of 
the  physical  quantity  on  the  plastic  and  rigid  region, 
respectively  (same  below).  Under  the  assumption  of  square  yield, 
the  orthogonal  flow  method  and  equations  (16)  and  (32)  in 
reference  [1]  are  used.  When  the  rigid  plastic  interface  is 
stationary,  if  the  stress  points  on  the  plastic  side  are  located 
on  AB  or  CD  of  the  yield  plane,  then  we  have 

[y]«0  and  [Q]=[M]=0  (8) 
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4 

If  the  stress  points  on  the  plestlc  side  are  located  on  AC  or  BD 
of  the  yield  plane,  then  we  have 

[*3-0  and  CQ]-tM]-0  (9) 


Figure  3 
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III.  Dynamic  Analysis 

Assuming  that  the  yield  plane  Is  square,  the  results  of 
rigid  plastic  limit  analysis  (taking  shear  effect  into  account) 
of  a  fixed-ended  beam  under  uniform  load  give  the  initial 
mechanistic  phase  diagram  as  shown  in  Figure  4(c).  When  u<0  and 
0<p<8,  the  beam  is  stationary.  When  0<&$8  and  p  |  t_Q=p(0)>& ,  the 
beam  slides  with  respect  to  the  support.  The  solution  is  easy  to 
find  and  it  is  independent  of  a.  When  0>8  and  p(0)>8,  there  is  a 
plastic  region  in  the  middle  of  the  beam.  The  stress  points  are 
located  on  side  AB(non-corner )  of  the  yield  plane.  The  region 
between  this  region  and  the  end  of  the  beam  is  the  rigid  region. 
When  the  beam  end  is  a  plastic  hinge,  and  stress  points  are  on 
side  CD,  if  stress  points  are  at  the  corners  C  and  D,  then  the 


but  it  tht  «nd  not  only  rotates  but  also  slides.  Otherwise,  the 
baeai  only  rotat«r~  In  the  following,  we  will  discuss  the 
situation  u(0)>8  and  ft>8.  Because  of  symmetry,  only  the  left 
half  of  the  beam  is  discussed. 
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Figure  4.  Initial  Mechanistic  Phase  Diagram  of  the  Beam  Under 
Three  Different  Considerations 

(a)  considering  bending  alone 

(b)  considering  shear  effect 

(c)  Timoshenko  beam 


Let  us  assume  that  y I  tatQ*y 1 *  T^e  b°un<*ary  conditions 
are;  M|x-0—  1,  Qlx,^»0.  Furthermore,  when  ya«y I x«o^° »  Qa=Qlx=0 
*1.  From  equations  C2),  (3),  (5)  and  boundary  conditions,  we  get 
the  following  solution:. 

In  the  rigid  region  (Osx^e): 
i  •  i,  +  d.*»  9  —  t,  +  <?«* 

PQ  -  PQ.  -  (2?  -  f.)x  + 

U  —  —  1  +  (p  —  •$.)*  — (2/m  —  +  7  i*** 

2  6 

In  the  plastic  region  (£<x<^): 


9  -  Ceh  X-^M2  +  2ft,  PQ- 

V  • 


i  *~U/2.  M  - 1 


(id 
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V* 


where  ^  represents  the  non-d leans  ional  value  of  the  angle  of 
rotation" in  the  rigid  region. 

(1)  When  y#«0  end  t<0,  by  aubatitutlng  equations  (10)  and 
(11)  into  the  discontinuity  condition  (6)  to  simultaneously  solve 
the  problem: 


where 


«.-[(•♦  4  f)  *  i==SA  -  VT  S«k •  ia 

IN  2  /  V  m  ^  J  B 

!  —  fit.-  -  */(A). 


From  ♦alt,0>0'  and  ?<0*  we  get 


where 


#40)- *!.«>••  j»>* 


A  -  2(«?  -  3fi  +  «•)/[«(«  +  <■)]• 

tfc  ~4/2  -  WTlMti  -  «)/(«?  -  «fl  -  12 «) 

V  • 


(12) 


(13) 
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(14) 
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Obviously,  When  u(0)<8,  y*0.  When  ,  Q  =1  and  y  >0. 

When  u«0 ,  5*0 . 

(2)  When  y  «0  and  5>0,  by  substituting  equations  (10)  and 
(11)  into  the  discontinuity  condition  (7),  we  get  the  solution: 


e-dZu-lt*- 

r-ls-itt-j^wo- --«>**]/ 

<15) 

<  *  [£  ”  ** + "2  ]/[(!.  **)  ^ + ,i)] 

$-{-*(*-  7/(«) + ^>(7  «*>  -  ct + ^*>  r.  **) 

-  «,*(l)  +  «XI»  •  (i'^)1}/[(C^r)C,,  +  ,,‘)] 

where  t=t-t0  ,  and  t,  is  the  initial  time  of  the  initial  state. 

*  “  *!»■•»  D—  j'z^fr  +  *({)  —  fr'(f)»  £  —  •  D/B  + 

[j  '*  ‘  [* ^  we> "  *«»'*  4+  y7."({)]/(j;  ** ), 

*** ~  *^7*  ” '^*'({)  ‘ d,t^7  J/B  6+  £ (,(o - «••({)){«. 

Notice  that  y|^=Q=0.  By  mathematical  deduction  we  can  prove 
that:  at  any  t#  we  have  v>0,  v’>0,  v"<0  and  v'x_^=0.  Hence,  m> 

0,  f(5)^0.  By  using  ia|5=%  =^al  ,,<;alr=t=o>0,  based  on  ^a*3  and  E>0. 

T=t=0  *  1 


can  obtain: 


0(0) > s,  0(0 - 7 j^r > *»»  **(*■) > * 


(16) 


19kM  »(0)<8,  fr»0.  When  iKt)»ut  ,  the  beam  reaches  Its  maximum 
deflection.  When  S(*)*u»  ,  t*0.  If  Qa»1,  we  get  *Ci ,  which 
is  a  contradiction  to  fc>0.  Therefore,  Qalt>o<1*  Furthermore,  ya 
-0  and  €*  C»  . 

(3)  Wien  ya»0  and  t»0,  we  will  first  prove  that  [$]»0  when  C 
■0  and  then  discuss  the  solution  in  this  situation. 

Let  us  assume  that  t«0  when  t^  £t<ta  and  when  t<ta  or  t<ta 
.  Note  that  A»2i»-ta  and  £|t-t^»£a  *  ^en  ^  <t<t*  ,  equations  (10) 
and  (11)  are  substituted  into  the  discontinuity  equation  (8)  to 
solve  C  and  $a ;  Then  we  get: 

Id!  «  {«.  -  *[«.  +  (•  -  y  fi)  *  fc=i/2/B.]/ 

(i7) 


where 


v  « 


If  t,>tt  and  ya*0,  we  get  |»0  when  t»  <t<t,  .  From  equations 
(12)  and  (15)  we  get  p||-0|lijji  ft^ijp  fc*0,  i.e.,  ***»•  I . 
0Qaobtained  from  equation  (8),  however,  is  only  a  function  of  u 
and  ?a  .  Hence,  &Qa“9Qa  t.t*6Qa|tMt  *  *n  t^le  two  cases  that  ta  > 
t,  ya>0  and  t^«t| ,  this  equation  is  obviously  valid. 

By  using  the  reverse  method  we  can  prove  that  if  y  >0  we 
have  t*0.  Hence,  when  t<t!  or  t>ta  ,  ya«0.  Its  solution  can  be 
expressed  by  equation  (12)  or  (15).'  When  £-0,  we  have 


■otlee  that 


Thtrcfort: 


[Ml  -  2  -  (fQ, -•*.);  +  -  0, 

0 


(,+(.-  ‘«j*k=jaw 


(18) 


By  substituting  equation  (18)  into  (17)  we  get  C4»)*0.  Thus,  we 
proved  that  discontinuity  condition  (6)  could  be  used  when  $»0. 
Hence,  when  ya-0  and  5*0,  the  solution  of  the  problem  is  the  same 
as  equatidn  (12).  Furthermore,  the  necessary  condition  to 
maintain  the  initial  state:  y(0)>8,  m<u,  and  w»0. 

(4)  When  ya>0,  Qa*1«  As  described  above,  C*0.  By  using 
discontinuity  condition  (6)  we  get 

I -I..  2(n-*)  (19) 

•  B,  B,  J » 

where  Bl  »C»,  sh  -*5A'ra-/a”  ch  Sj-V'TST.  Note  that  T=t-t,  and  t,  is 
the  starting  time  of  the  initial  state.  From  I  T3_o>® *  ^aU=0=® 
and  7al t>0>0*  we  «et: 

*U>*. 


When  S»i»x  ,  ya*0,  ya$0.  It  was  pointed  out  earlier  that  when  ?a=0 
we  have  CH*  .  Afterwards,  it  enters  the  state  of  yfl=0  and  e>0. 

The  above  discussion  is  applicable  to  problems  with 
uniformly  distributed  initial  velocity  as  long  as  we  make  y  (())-*• 

**  l  t>o  *°  and  Jo2**dt“ylt«o- 


IV.  Results  and  Discussion 
Table  1.  External  Total  Impulse  Load  (10  lb  s) 


«i  (w*0 


•*** 


Ca_ 

L_, 

0.1X19* 

0.1X10* 

0.3X10* 

0.2X10"* 

0.2x10"* 

0.2X10"* 

1.90071 

2.04200 

1.42113 

IU. 
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load  curve 

maximum  load  (Ib/in) 

duration  during  which  load  is  applied  (s) 
maximum  deflection  in  the  middle  (in) 
initial  velocity  problem 
initial  velocity  is:  0.2732(in/s) . 
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Figure  5.  Middle  Beam  Deflection  vs.  Duration  of  Loading  When 
Total  Impulse  and  Loading  Mode  Are  Identical 

1.  initial  velocity  problem  y/._n=0.2732  x  10“^in/ 


Figure  6.  Maximum  Deflection  Middle  and  End  of  Beam  vs.  a 

1.  maximum  limit  of  deflection  in  the  middle  of 
Timoshenko  beam  when  a~0 

2.  maximum  deflection  at  the  middle  of  the  beam 
according  to  classical  solution 

3.  maximum  middle  deflection  of  Timoshenko  beam. 

4.  maximum  end  point  displacement  of  Timoshenko 
beam  (beam  parameters  are  same  as  those  in 
Figure  12) 


1.  From  the  results  tabulated  in  Table  1  we  can  see  that 
the  maximum  deflection  at  the  middle  of  the  beam  is  different 
when  the  total  external  impulse  load  is  the  same  but  the  loading 
mode  is  different.  Figure  5  shows  the  middle  beam  deflection  vs 
loading  time  curve  when  the  total  external  impulse  and  the 
loading  mode  are  identical.  We  can  see  that  when  the  loading 
time  approaches  zero,  the  maximum  deflection  at  the  middle  point 
of  the  beam  approaches  the  corresponding  value  of  the  initial 
velocity  problem. 

2.  Rotational  inertia  and  shear  effect  have  an  obvious 


effect  on  the  mechanism  of  the  motion  of  the  beam.  Figure  4 
shows  the  initial  mechanistic  phase  diagrams  of  the  fixed-ended 


beam  under  the  different  considerations  to  explain  this  effect. 

In  Figure  4(c),  the  curve  |i(0)*f  (a,$)  separating  the  rotational 
mechanism  from  the  rotational  and  sliding  mechanism  is  given  in 
the  following  equation  (they  are  obtained  by  changing  jix  to  p(0) 
in  equation  (14)). 

*0)  -  2(«!  -  +  6m)/[gKpt  +  ««)] 

lb  ~1A  -  */^.(>f.  -  4)/(«i  -  “  12.) 

V  « 

3.  Effect  of  Rotational  Inertia  on  Maximum  Deflection 
In  reference  [6],  Jones  showed  Figure  7  to  explain  the 
effect  of  rotational  inertia  on  the  maximum  deflection  of  the 
beam  in  initial  velocity  problems.  In  this  paper,  after 
considering  the  load,  the  results  are  shown  in  Figures  6  and  7. 
From  Figure  6  we  can  see  that  under  a  specific  load,  without 
changing  other  beam  parameters,  the  larger  the  value  of  a  is,  the 
smaller  the  mid  point  deflection  is.  When  a  increases  to  a 
specific  value,  shear  slide  begins  at  the  end  point  of  the  beam. 
When  o  approaches  zero,  the  maximum  deflection  of  the  middle  of 
the  beam  approaches  a  limiting  value.  However,  this  value  is 
higher  than  the  maximum  deflection  given  by  the  classical 
solution  in  which  only  the  effect  of  bending  is  considered.  This 
situation  can  be  explained  as  follows: 


firm  (iW<) 


Figure  7.  Rigid  Plastic  Interface  Plane  Positions  i  and  Mid 
Beam  Velocity  vs.  Time  with  Different  Rotational 
Inertia 

Beam  parameters:  M,  =0.191  x  10  Ibin;  Q0 =0.880  x 
10*  lb,  m=0.366  x  10_*Ibsa/in,  l=5in,  6=0.231  x  10*  , 

.J 

(when  o=0.8329  x  10  ,  it  is  equivalent  to  a  square 

cross-section  beam  with  uniform  mass  distribution, 
o0  =30.5  x  10*  lb/ ina  ,  t0=o0|vT  ,  p=0.732  x  10'*  lbs4 /in* 
h=0.5  in  and_b=1in,  where  o#  ,  t„  ,  and  p  are  the  unit 
tensile  yield  stress,  shear  yield  stress  and  density 
of  the  beam  material,  and  h  and  b  are  the  height 
and  width  of  the  beam  cross-section). 

1.  classical  solution 

2.  classical  solution 


From  Figures  7(b)  and  (c)  we  can  see  that  when  a  is  very 

small,  £,  the  time  to  reach  maximum  deflection  t^  and  the  mid 

point  velocity  y|  ,  when  t<t_  given  by  this  solution  and  those 

x=;?  s 

given  by  the  classical  solution  are  very  close.  Here,  t 

S 

represents  £=%  in  the  classical  solution.  However,  when  t>t2  , 
the  values  of  y|Xs_^  given  by  the  two  solutions  are  quite 
different. 

This  solution  gives  / 51 

“  Ct"  +  o  >  O,  \  «-*o. 

Let  a-0,  we  get 

The  classical  solution  gives  y/x=^=*;DT2iidt-24(t-ts)  , 
and  (t>t  ).  Here,  t  is  the  duration  over  which  the  external  load 
is  applied.  It  is  obvious  that  the  former  is  a  uniform  motion 
and  the  latter  is  a  uniform  deceleration  motion.  These  two 
equations  are  integrated  to  obtain  the  maximum  deflections  of  the 
two  solutions  at  the  middle  of  the  beam.  In  this  example,  they 
differ  by  36.3%. 

This  work  was  completed  under  the  guidance  of  Professor  Wang 
Ren.  Our  teacher  Huang  Zhuping  also  assisted.  The  authors  wish 
to  express  their  gratitude. 
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AN  ANALYTICAL  SOLUTION  OF  DYNAMIC  RESPONSE 
FOR  THE  RIGID  PERFECTLY  PLASTIC  TIMOSHENKO  BEAM 

Jin  Qhanlin 

(Smart*  Imtkau  of  Mtchsmcai  mid  EUcthcd  Ttthmoiogy  Ministry  of  hUchitu-B*iUi*i  Industry) 


Abstract 

Pmfoos  studies  on  Timoshenko  beam  was  only  dealt  with  for  impulmve  Telocity 
but  not  fur  iaereaing  load.  An  analytical  solution  of  dynamic  response  for  the  fixed- 
ended  Tbatshanko  beam  subjected  to  uniformly  distributed  dynamic  load  is  (inn  he¬ 
rein  by  see  of  the  discontinuity  conditions  for  moving  rigid-plastic  The  solu¬ 

tion  is  valid  for  arbitrarily  time-varying  but  non-reversing  load.  Finally  the  influence 
of  rotatory  inertia  on  the  dynamic  plastic  response  of  beams  is  discussed. 


General  Variational  Theorem  for  Structural  Plastic  Buckling  ' 
Analysis  Using  Deformation  Theory 
Li  Guochen 

(Institute  of  Mechanics,  Academia  Sinica) 

Abstract 

This  paper  gives  a  class  of  general  variational  theorem  to 
analyze  structural  plastic  buckling  using  deformation  theory, 
which  explains  the  essential  significance  of  potential  energy. 

In  the  form  of  general  variation,  we  proved  that  there  is  no 
unloading  during  plastic  buckling.  Finally,  it  was  applied  to 
examples  'in  the  analysis  of  reinforced  plate  and  shell. 

Introduction 

As  we  know,  in  the  structural  analysis  of  plastic  buckling, 

it  is  more  appropriate  to  use  the  increment  theory  of  Prandtl- 

Reuss  to  describe  'the  stress-strain  relation.  The  results  are 

even  more  close  to  the  experimental  data. 

Regardless  of  whether  it  is  elastic  buckling  or  plastic 

buckling,  the  analysis  must  be  carried  out  in  two  parts.  One  is 

to  find  the  basic  solution  prior  to  buckling.  Second  is  to  see 

whether  there  is  another  path  at  any  point  on  the  basic  route. 

Ill 

When  conducting  the  buckling  analysis,  Hill  J  used  a 

"comparative  elastic  solid"  model,  i.e.,  there  is  no  unloading  at 

[21 

the  instance  of  buckling.  Hutchinson  analyzed  the  second 

order  generalized  function  consisting  of  displacement  variables 

til 

based  on  the  uniqueness  of  Hill's  solution  to  further  verify 


this  problem.  Bushnell  J  used  a  difference  method  to  calculate 
many  structural  buckling  problems  and  he  claimed  that  the  no 
unloading  phenomenon  is  a  condition  for  "consistent  loading". 

The  history  of  this  condition  was  also  introduced  in  reference 
[3] . 

The  computation  of  plastic  buckling  is  far  more  complicated 
than  that  of  elastic  buckling.  In  order  to  simplify  the  problem 
to  the  extent  possible, the  variational  principle  is  a  powerful  tool. 

[A  ] 

Referring  to  the  generalized  variational  theorem  under  elastic 
conditions  developed  by  Reissner,  this  paper  gives  a  general 
variational  theorem  in  the  analysis  of  structural  buckling  using 
deformation  theory.  The  physical  significance  of  the  generalized 
function  and  its  second  order  variation  involved  is  explained. 

The  basis  of  "consistent  loading"  used  in  this  general 
variational  theorem  was  further  proved.  Finally,  the  superiority 
of  this  general  variational  theorem  was  demonstrated  in 
simplifying  the  basic  equations  of  plate  and  shell  structures. 


'*  V  V 
•  >  * 


•  V  V 
»  "  •  *  « 

«  *  •  * « 


I.  A  Class  of  General  Variational  Theorem 
Let  us  assume  that  the  generalized  potential  energy  is 


v  t\ 
V  v 


n  -  ^  -  j  T.mtdfr 


where  is  the  stress  tensor,  is  the  strain  tensor  which  is 
related  to  displacement  according  to  the  following  known  relation 

««  —  4-  (■«  + 

2  2 


i  V.  ■-  *  * 


•  ■*  .'•* 
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Here,  the  comma  in  front  of  a  subscript  represents  the 
differentiation  with  respect  to  that  subscript.  is  the  known 
externally  applied  stress  on  the  surface  s^,. 
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Here 


Vr(3/2>(>*«,  »u 


When  oc<oT 

4**0 

E,  v,  Eg  and  or  are  the  elastic  modulus,  Poisson  coefficient, 
cutting  modulus  in  unit  tensile  and  yield  stress,  respectively. 
The  first  variation  of  TT  is 


"  - [.»».„  +  (.« -  *  -  J,  ***« 


(2) 


where 

8ta  -  —  (*•..<  +  «■;.,)  +  -7  +  ■*..*»«) 

2  2 

Using  divergence  theorem  we  can  derive  that 

+  (»»«•<.«)#  i}8u/dw 

+  U*.!  —  [-jCO  +  *V.»  -  1„  ||  so-i,* 

+  i  U*.v  +  (*,>*„«)]"#  —  ft}auidsT 

J'T 

+  I  {»«  +  (»«•„*)}»»  Midi,  —  0 

*  *m 


(3) 


where  Oj  is  the  normal  component  of  the  reference  coordinate 
prior  to  deformation  on  the  boundary.  It  is  not  difficult  to  see 
that  the  Euler  equations  in  (3)  corresponds  to  the  well  known 
equilibrium  aquation,  constitutive  relation  in  deformation  theory 
and  boundary  condition. 

Equation  (1)  is  similar  to,  but  different  from,  those  of 

[A] 

Reissner's  in  an  elastic  system.  If  we  continue  to  proceed 
with  a  second  order  variation  and  notice  that  u^  and  are 
Independent  variables,  «* u^«a* j»0.  Based  on  these,  we  get 

*‘B~  Lb1*** +'"*‘-+(***_5&*'*H*  (4) 

where 

•*««  ■»  *■*«*■*#. 

In  the  following,  a  new  variation  «*  (refer  to  reference 
C8])  is  made  with  respect  to  6o^j  and  6u^  in  the  second  order 
generalized  function  Q  to  obtain  relevant  equations  and  boundary 
conditions  during  buckling.  Noticing  the  "consistent  loading" 
condition,  the  stress,  displacement  and  plastic  modulus  remain 
unchanged  in  the  variation  a*  process.  The  "consistent  loading" 
condition  and  its  verification  in  general  variation  will  be 
discussed  in  the  following. 

By  using  the  divergence  theorem,  we  can  derive  the  following 
from  equation  (4) : 


•  “2  4"  (toft*, 4),  j  +  (ffnSmt.0,  i}8*{8»t)dp 

♦  2  ^  {#•*  -  [i  ((1  +  p)*r«  -  *«***,,) 

+  +  J  **«  ]}  8*($a,ddw 

+  2  {(r«  +  ■  ® 


E,fc  la  the  tangential  modulus  during  unit  tensile  stress.  The 
Euler  equations  In  (5)  correspond  to  the  known  equilibrium 
equation  during  buckling,  incremental  constitutive  equation  of 
deformation  theory  and  boundary  condition  during  buckling. 

Substituting  the  relation  between  6*^,  6**^  and  du^  into 
equation  (4)  we  can  also  get 

£>  —  —  /  +  (*»**«<.«)»  ' 

+  {*««  —  [-J  ((l  +  »)»*«  —  **»«*«)  ^ 

+  <0  +  —  8s»  | J  8a,,  dt 

4*  |  {9fa  +  80n*i,i  +  alk8mi^)n)8uiit 

Comparing  equation  (6)  to  (5)  we  can  see  that  Q  and  6*Q  are 
similar  in  form^This  is  the  same  in  the  elastic  case. 
Therefore,  when  fi*Q»0,  Q»0  and  vice  versa. 

Thus,  we  proved  that  the  first  order  variation  of  the 
generalized  function  TT  in  equation  (1)  can  be  used  to  derive 
various  basic  equations  in  the  basic  path  to  solve  the  problem 


A 


prior  to  buckling.  By  carrying  out  a  new  variation  6*  with 
respect  to  its  second  order  variation  can  result  in  relevant 
equations  during  buckling. 


II.  Physical  Significance  of  Generalized  Functions  TT  and  Q 
Now,  let  us  decompose  the  first  term  in  the  integral  in 
equation  (1) 


.H.  -  p'V...,)  -  + pv*. 


(7) 


where  the  function  to  be  integrated  e^j  can  be  expressed  in  the 
form  of  stress  and  can  be  expressed  in  the  form  of  strain. 
The  second  term  of  equation  Cl)  is 


V 


<*«>  QV 

!• 


A?* 


(8) 


When  ejj-aV/aoj^  is  valid,  we  get  the  following  by  substituting 
equation  (7)  into  it 

Stt  n> 

,  *lid9H  '  (9) 


The  right  term  represents  the  strain  energy  per  unit  volume. 

Thus,  TT  in  equation  (1)  is  naturally  equivalent  to  the  definition 
of  potential  energy.  For  this  reason,  at  the  beginning  of  this 
paper  we  called  it  the  generalized  potential  energy. 

In  order  to  determine  stability,  the  classical  definition  / 5 1 5 
proposed  by  Dirichlet  and  Kelvin^^  is  that  the  body  is  in  a 
state  if  the  total  amplitude  of  the  added  displacement  caused  by 
the  deflection  at  any  instance  is  arbitrarily  small  when  the 
deflection  itself  is  arbitrarily  small.  On  the  contrary,  if  an 


arbitrarily  snail  deflection  would  lead  to  a  finite  amplitude 
change,  then  it  is  unstable.  Obviously,  the  necessary  and 
sufficient  condition  is  that  the  amount  of  internal  energy  stored 
or  consumed  when  an  infinitesimal  displacement  is  added  to  the 
equilibrium  position  must  be  larger  than  the  work  done  by  an 
external  force.  In  the  1930's,  Trefftz^  and  Kappus^®^  also 
applied  this  principle  to  derive  the  elastic  buckling  equation 
and  stability  determination  criterion  in  the  form  of  energy.  In 
conjunction  to  the  topic  of  this  paper,  for  an  arbitrarily  small 
deflection  (AOjj,  Au^),  if  Air>0,  it  is  stable,  if  Aw<  0,  it  is 
unstable. 

-  From  Taylor  series  expansion  we  know  that 


n  +  ah  -  [(*„  +  *r„)  +  ut,  +  ^  «*•*) 

-  ( r + ^  + « * 

—  (  ?,(«,  + 

*• T 


Here,  the  external  stress  T^  does  not  vary  with  deflection. 
Therefore,  it  is  a  case  of  "fixed  load".  By  expanding  the  above 
equation,  we  get 


H  +  AD  -  t <r»*H  ~  FI*  - 

+ + (•«  ~  i^)  H  *  ~  Lr  7'Suiitr 

+  {i  Jr[  +  Mt#  +  (»'« -  a^T *a*J) H dV 


(10) 


4. ...  -  n  +  sn  +  -r-fn- 


Ncauit  Air  is  already  equal  to  zero,  therefore,  whether  it  is 
stable  depends  on  A*ir(«Q).  if  Q>0,  it  is  stable.  If  Q<0,  it  is 
unstable.  Q«0  is  the  stability  limit.  Above  the  stability  limit 
or  on  buckling  points,  the  criterion  for  stability  is  determined 
by  higher  order  of  variations  of  ir ,  such  as  Asir....  This  says 
that  Q»0  may  be  stable  or  may  be  unstable.  Therefore,  Q>0  is 
only  the  sufficient  criterion,  but  not  the  necessary  criterion, 
for  stability.  After  realizing  the  significance  of  Q-0  and 
knowing  Q»0  corresponds  to  A*Q=0,  the  reason  why  the  buckling 
equation  can  be  derived  from  A*Q*0  is  clear. 

In  addition,  the  following  explanation  is  given  to  describe 
the  -relation  between  the  condition  Q=0  and  the  loss  of  uniqueness 
of  the  solution.  If  the  uniqueness  of  the  solution  on  the 
loading  path  is  lost,  then  there  must  be  other  possible 
increments  ao^  and  AUj^  which  satisfy: 

Aa,H  +  i  +  ;  -  0 


A 


—  4r  [(1  +  v)AOi,  —  vA + 


(11) 

(12) 


AT,  —  (Affii  +  Aa,k.MM  +  ffnA»,jl')»i  ”  o 

Am,  —  0 


on  s. 


or 


on  s 


u 


(13) 


If  we  choose 

Affn  ■  Am,  ■  Sm, 

We  can  see  that  the  Euler  equations  in  (11)— (13)  are  similar  to 
those  in  (5).  This  means  that  the  uniqueness  of  the  solution  is 
lost  during  buckling  and  vice  versa.  Therefore,  the  buckling 
point  must  be  related  to  the  loss  of  uniqueness  of  the  solution. 
However,  stability  is  not  related  to  uniqueness,  because  buckling 


Is  not  equal  to  unstable 


In  simplicity,  buckling  and  loss  of  uniqueness  of  the-  / 5 

solution  correspond  to  Q»0  or  d*Q*0.  Whether  it  is  stable 
depends  on  Q>0  or  Q<0.  At  Q-0,  it  depends  on  the  sign  of  higher 
order  variation  terms  in  equation  (10). 

III.  Investigation  of  "Consistent  Loading"  Condition 
Hutchinson ' sI4argument  can  be  extended  to  generalized 
variation  in  the  following  steps. 

Under  a  "fixed  load",  the  amplitude  variation  of  the  force 
(or  displacement)  load  is  proportional  to  a  parameter.  Its 
distribution  is  Independent  of  this  parameter.  It  was  explained 
earlier  that  the  necessary  and  sufficient  condition  for  stability 
is  Q>0.  Now,  let  us  re-write  equation  (4)  into 

Q  —  +  ««**•*.<)  + 

—  [-j  (O  +  »)•*«  —  «*»»«•«)  (14) 

+  «  +  7*")]*"'} 4V 

where 

1  plastic  loading 
a  »^0  elastic  loading 

0  plastic  loading  followed  by  elastic  unloading 
Let  us  take  a  solid  body  for  comparison,  i.e.,  an  elastic 
comparative  solid  body.  Its  second  order  generalized  function  is 
defined  as 

Qc*right  hand  terms  in  (14),  but  s  *  p  plastic  loading 

(0  elastic  loading 


It  Is  obvious  that  there  is  no  elastic  unloading  zone  In  this 


case.  To  this  end,  an  infinitesimal  deflection  (60^,  6  u^)  Is 

added  to  a  known  based  state  C j  ,  u^).  The  difference  between  Q 

and  0  is 
c 


where 


(15) 


(16) 


The  subscript  c  represents  the  possible  presence  of  an  elastic 
unloading  zone.  By  substituting  the  expressions  for  4  and  64 
into  (16)  we  get 


because 


'■ ~  *L  Id,  ~  + f  (7.  ~ 


therefore  Ic^0 ,  Q>QC . 

This  is  to  say  that  Q=QC  if  a  certain  deflection  is  added  to 
the  entire  body  area  V  so  that  the  plastic  region  remains  plastic 
and  elastic  region  remains  elastic.  Otherwise,  Q>QC  if  the 
additional  deflection  leads  to  the  presence  of  an  elastic 
unloading  zone  in  V  .  Therefore,  it  is  only  necessary  to  find 
the  corresponding  characteristic  function  (60^,  6u^)  in  Qc  in 
order  to  determine  the  minimum  critical  loading  parameters 
because  it  is  not  possible  to  have  even  smaller  critical  loading 
parameters  and  other  characteristic  function  distribution  to  make 
Q«0  while  keeping  Qc>0.  ’  Thus,  we  proved  that  the  "consistent 
loading"  condition  is  still  valid  in  the  form  of  generalized 


variation. 


IV.  Application  in  Analysis  of  Plate  Shell  Structure  / 5 1 7 
In  a  thin  wall  plate  shell,  let  us  assume  that  the 
distributions  of  strain,  stress,  plane  displacement  and  plastic 
parameter  along  the  x^-direction  are  approximately: 


-  [«S.  Njkt  8t%  8NJk,  uT,  **'] 

•+  (*.#,  12MJV,  8K*,  128MJV,  2^/iU,  (17> 

where  a,  6=1  or  2  and  h  is  the  wall  thickness.  By  substituting 
the  above  into  equation  (1)  and  using  the  strain-stress  relation 

[O] 

in  a  plate  shell  ,  we  can  get  the  following  by  integrating 
along  and  applying  the  divergence  theorem: 


an  -  \'  {tAU.]*.? 

+  [  ' “  Ek  (  + 1 c V**)] 

+  [a-  -  ^  (■ j  C*irtNrt  +  C^Mh)]  *  . 

+  (  [•••]^Lr+(  I  •  *  *  —  0 

)tr  Jim 


(18) 


where  w^^  is  the  initial  deflection  and  bog  is  the  curvature 
tensor 

CSftrf  -  (1  +  *)  [y  (a„8M  +  a«a, T)  - 

+  <A"’  +  «-*#,)  -  '  (19 

CSSh  -  [*7  (A.r*«  +  -  J  Wrl] 


116 


•?.  • . 


Similar  to  the  algorithm  used  in  the  first  section,  we  apply 
a  new  variation  5*  to  n  for  a  second  time.  Then,  we  have 

ana’ll)  -  fQ 

-  -2  (  {[aA^la*^) 

•  • 

+  18M**  +  8N.+*  +  8NJ.»  +  wu>).  w 

+  *l8\8m)Us 

+  2  j  {[  8t*i  -  ±  D^t8Mrt )]  S*(»lVw) 

+  [a*„  -  (|  BSL8N+  +  OiW"r.)]  **('*->}  A  (2Q) 

+  ••• 

where 

DS X*  —  (1  +  »)  |  y  (8^Bm  +  8*8*,}  —  ^  +  ~  8-*8r*  j 

+  (8„8m  +  8«8»)  -  i-a^ri]  +  *Sr. 

^  [-J (8„am  +  a^)  - y  a.**]  +  *2U  (21 ) 

/  5 1 8 

**"  "  *4  (£,  “  f )  Cwfr#  *  ^  +  7 
£•*  "■  ~  [  y  +  arta^)  —  y  a)l#aT«J  v* 

From  the  above  it  is  obvious  that  the  third  and  fourth  Euler 
equations  in  (18)  and  (20)  correspond  to  the  constitutive 
equations  before  and  after  buckling.  From  (19)  and  (21)  we  can 
see  that  the  rigidity  matrix  correlating  generalized  strain  to 
generalized  stress  is  symmetric.  However,  when  another  method 
was  used  to  derive  it  in  reference  [10]  the  favorable  condition 
of  matrix  symmetry  was  not  realized. 


In  the  above,  the  plate  shell  belongs  to  the  case  of 
symmetric  upper  and  lower  cross-section  relative  to  a  neutral 


plane.  For  a  non-symmetric  cross-section  with  reinforcing  ribs, 
an  appropriatley  chosen  generalized  force  can  be  used  to 
simplify  the  equations. 

Let  us  assume  that  there  is  a  rib  along  the  x-direction 
(neglecting  twisting  resistant  rigidity).  Its  cross-sectional 
width  is  b(z),  i.e.,  it  can  vary  along  the  z-direction.  The 
stress  and  strain  distribution  on  the  cross-section  are 


(22) 


where  A  and  I  are  the  cross-sectional  area  and  the  moment  of 
XX. 

inertia  relative  to  the  neutral  plane  of  the  ribbed  cross- 
section,  respectively. 

After  substituting  them  into  (1)  and  through  certain 
operations  we  get 


(23) 


N.*,* 
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where  H  is  the  height  of  the  rib  bar  itself.  Again,  we  have 
.  A .  -  j  Mi,  J  -  \m*,  I  -  0  -  $/«*'* 


Let  us  assume  that  the  spacing  between  the  origin  of  the  z 
coordinate  and  the  neutral  plane  of  the  rib  bar  is  e  (usually  the 
original  can  be  chosen  at  the  center  plane  of  the  plate  shell), 


hence 


j  —  A,t,  1  —  /,  +  A,f 


From  equation  (23),  in  order  to  simplify  the  equation  to  the 
extent  possible,  especially  in  the  equilibrium  equation,  we  can 
define  a  new  generalized  force  as: 


Thus,  we  can  solve 


Substituting  back  into  (22),  we  get  / 5 1 9 


Based  on  the  above  we  know  that  although  the  relation 
between  stress  and  generalized  force  is  complicated  by  choosing  N 
and  M  as  the  generalized  forces  (such  as  (24)  in  comparison  to 
(22)),  yet  the  equilibrium  equation  can  be  simplified. 


V.  Conclusions 


VV 


The  variational  theorem  introduced  to  analyze  the  structural 
plastic  buckling  based  on  deformation  theory  not  only  can  be 
applied  to  mixed  boundary  problems  but  also  provides  a  way  to 
simplify  the  equations  to  be  solved. 
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A  GENERAL  VARIATIONAL  THEOREM  FOR  THE 
'STRUCTURAL  PLASTIC  BUCKLING  ANALYSIS 
USING  THE  DEFORMATION  THEORY 

Li  Guochen 

(fmuitmtt  tf  Mttkamitt,  KctdtmU  Simtt) 


Following  the  farm  given  by  Reismcr  in  19S0  for  elude  analysis,  a  general  variational 
functional  in  plasticity  is  prescribed  u  IT  in  eq.  (I).  By  setting  its  first  variation  311  due  to 
the  variation  of  (<fy,  *<)  to  zero,  the  Euler  equations  derived  in  (3)  are  proved  to  he  the 
equilibrium  equations,  a  deformation  type  of  stress-strain  relations  and  boundary  conditions  for  the 
prebuckling  fundamental  path  solution.  As  Kappus  bad  done  in  1939,  a  new  variation  5*  can 
be  imposed  on  (*r, *,*»,).  Let  Q  —  d*U,  then  from  B*Q  of  eq.  (5),  the  basic  incremental 
equations  are  derived  for  die  evaluation  of  critical  loading  and  its  corresponding  budding  patt¬ 
ern.  Comparing S)  and  (6)  it  can  be  seen  that  whenever  8*Q  equals  zero  the  same  is  Q 
or  vice  vena. 

Form  eqs.  (7)— (9)  it  is  dtown  that  eq.  (1)  is  essentially  equivalent  to  the  potential 
energy.  According  to  the  definition  of  Dirichlet  and  Kelvin,  stability  depends  on  whether  AU 
is  positive  or  negative.  Using  the  Ttylor  series  expansion  eq.  (10)  brings  out  that  (a)  if  Q 
>  0,  stable  (b)  if  Q  <  0,  unstable  (c)  when  Q  *•  0,  buckling  occurs,  it  is  the  limit  of 
stability,  the  stability  at  this  point  .relies  on  the  sign  of  the  higher  variation  term,  e.  g.  3*11  •  •  • . 
On  the  other  hand,  when  the  uniqueness  of  die  solution  fails,  then  the  possible  incremental  para 
(AoVi  and  A*,-)  should  satisfy  eqs.  (11) — (13),  which  are  similar  to  the  Euler  equations  in 
(5). 

A  comparison  solid  which  has  no  unloading  condition  within  the  plastic  region  at  the 
moment  of  buckling  is  introduced  to  solve  the  buckling  problem.  Application  of  the  above 
theorem  in  the  plete  end  shell  problems  is  ezemplified. 


